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Abstract 



The software tool GRworkbench is an ongoing project in visual, numeri- 
cal General Relativity at The Australian National University. This year, 
GRworkbench has been significantly extended to facilitate numerical exper- 
imentation. The numerical differential geometric engine has been rewritten 
using functional programming techniques, enabling fundamental concepts 
to be directly represented as variables in the C++ code of GRworkbench. 
Sophisticated general numerical methods have replaced simpler specialised 
algorithms. Various tools for numerical experimentation have been imple- 
mented, allowing for the simulation of complex physical situations. 

A recent claim, that the mass of the Milky Way can be measured using 
a small interferometer located on the surface of the Earth, has been inves- 
tigated, and found to be an artifact of the approximations employed in the 
analysis. This difficulty is symptomatic of the limitations of traditional pen- 
and-paper analysis in General Relativity, which was the motivation behind 
the original development of GRworkbench. The physical situation pertaining 
to the claim has been modelled in a numerical experiment in GRworkbench, 
without the necessity of making any simplifying assumptions, and an accu- 
rate estimate of the efi'ect has been obtained. 
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Chapter 1 
Introduction 



GRworkbench is a numerical, visual tool for exploring analytic space-times in 
General Relativity. This year, the numerical differential geometric engine of 
GRworkbench has been rewritten using functional programming techniques, 
with the objective of creating a general platform in which complex physical 
situations can be simulated in numerical experiments. New tools for mod- 
elling physical systems were implemented within the functional framework. A 
recently proposed experiment, to determine the mass of the Milky Way, was 
analysed, and then investigated numerically in GRworkbench, demonstrating 
the applicability of the new techniques for numerical experimentation. 



1.1 Summary of thesis 

GRworkbench arose from work in visual numerical relativity by S. M. Scott, 
B. J. K. Evans, and A. C. Searle, at The Australian National University. Most 
recently, A. C. Searle implemented a numerical differential geometric engine, 
and improved 3-D visualisation [H]. The efficacy of the differential geo- 
metric engine, and the utility of GRworkbench as an intuitive visualisation 
tool, has been demonstrated [THE]. Chapter [2] presents an overview of the 
GRworkbench project. 

In order to facilitate the creation of a general system for numerical exper- 
imentation in analytic space-times, the numerical and differential geometric 
aspects of GRworkbench have, this year, been rewritten using functional pro- 
gramming techniques. Functional programming allows functions, like normal 
data, to be stored in program variables and manipulated by other functions. 
Important concepts in differential geometry, which are naturally thought of 
as functions, can thus be directly represented in the C-I-+ code of GRwork- 
bench. The functional programming methods employed in GRworkbench are 
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introduced in Chapter El 

Some of the numerical methods previously employed by GRworkbench 
were found to be too inflexible or inaccurate to be applied to potentially 
complex and computationally intensive numerical experiments. Sophisti- 
cated new algorithms have been implemented this year for key numerical 
operations including differentiation, integration, and minimisation; these op- 
erations act directly on functions, using the new functional framework of 
GRworkbench. A general notion of approximate equality permits the numer- 
ical methods to be implemented in a consistent and elegant way. Numerical 
methods are the topic of Chapter HI 

Appendix |X] lists the C++ code, written by the author, for the new 
numerical algorithms discussed in Chapter |H 

The differential geometric engine of GRworkbench, which relies on nu- 
merical methods for operations such as the transformation of tangent vector 
components between coordinate systems, has been rewritten within the func- 
tional framework, to interact cleanly with the numerical engine of GRwork- 
bench. Abstract notions such as points and tangent vectors are represented by 
C++ classes, which provide routines to obtain the coordinates of the objects 
in any coordinate system. The functional numerical differential geometric 
engine is described in Chapter [5l 

Physical situations in numerical experiments are modelled in terms of 
important objects in differential geometry, particularly points, tangent vec- 
tors, and geodesies. The key operation of geodesic tracing from initial data 
has been re-implemented using the new functional numerical engine. New 
methods for locating geodesies that are implicitly defined by boundary condi- 
tions have been developed using the function minimisation algorithms. These 
tools facilitating numerical experimentation in GRworkbench are the topic of 
Chapter [6l Appendix [B] lists the C++ code, written by the author, for a 
numerical experiment described in Chapter [HI 

An analysis of a recent claim by Karim et al. [8j, that the mass of the 
Milky Way can be determined using a small interferometer located on the 
surface of the Earth, is presented in Chapter [71 Properties of the interferom- 
eter model employed in the calculation of Karim et al. are investigated. The 
claimed size of the effect is found to be due to the coordinate-dependent def- 
inition of the interferometer employed, and not to the effects of space-time 
curvature. A more physically motivated interferometer model ('geodesic- 
defined interferometer') is proposed, and its properties are investigated. 

The interferometer model of Karim et al. and the geodesic-defined inter- 
ferometer were each simulated in GRworkbench. The results of these numer- 
ical experiments are presented in Chapter [HI The analysis by Karim et al. of 
their proposed interferometer was found to be in agreement with the results 
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of the GRworkbench simulations of that interferometer. The behaviour of the 
geodesic-defined interferometer was characterised, and used to obtain a new, 
more accurate, estimate on the size of the effect described in [8]. The effect 
was found to be too small to detect with an interferometer on Earth. We 
conclude that the proposed experiment, to measure the mass of the Milky 
Way using an interferometer located on Earth, is not currently technically 
feasible. 



Chapter 2 
GRworkbench 



GRworkbench is a software tool for visualising numerical operations on ana- 
lytically defined space-times in General Relativity. It has arisen from work in 
visual numerical relativity by S. M. Scott, B. J. K. Evans, and, most recently, 
A. C. Searle. In this chapter we give an overview of the motivation behind, 
and history of, the GRworkbench project. 

2.1 Motivation 

Analytic results in General Relativity are, in general, difficult to obtain. 
Exact solutions of the Einstein field equation are rare, and some physically 
important exact solutions are sufficiently complicated to be difficult to work 
with algebraically. It is usually necessary to make approximations if algebraic 
results are desired; this is exemplified by the claim analysed in Chapter [3 

Computational methods have been applied to the solution of the Ein- 
stein field equation for various boundary conditions, most famously to the 
currently unsolved problem of two in-spiralling black holes. Symbolic al- 
gebra software such as Mathematica, as well as specialised packages, such 
as GRTensorll and Sheep, are used to manipulate the tensor equations of 
General Relativity. 

Computational methods have also been used to explore the physical prop- 
erties of analytic solutions to the Einstein equation, through numerical opera- 
tions such as geodesic tracing. Traditionally, such simulations were performed 
using specialised codes as required. 

Visualisation in General Relativity is intrinsically difficult because space- 
times are 4-dimensional and curved, whereas computer monitors (and most 
other visualisation devices) are 2-dimensional and fiat. Traditionally, visu- 
alisation is performed by choosing a coordinate system, suppressing 1 coor- 
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dinate, and plotting the remaining three coordinates via a projection from 3 
dimensions to 2 dimensions. 

The goal of the GRworkbench project is to create a visual software tool for 
numerical General Relativity, in which a point-and-click interface encourages 
the user to explore freely in a space-time. Such a tool would, for the first 
time, allow experimental techniques to be applied to problems in General 
Relativity in an intuitive, visual environment. 

2.2 GRworkbench 

Working with S. M. Scott and B. J. K. Evans, A. C. Searle implemented a 
new version of GRworkbench in 1999 [15j. It featured an imbedded platform- 
independent GUI (Graphical User Interface), a novel numerical differential 
geometric engine, and a flexible visualisation system, as well as being easy 
to extend with additional space-time definitions. 

The differential geometric engine of GRworkbench allowed for abstract 
objects, such as points and tangent vectors, to have multiple numerical rep- 
resentations, corresponding to different coordinate charts. GRworkbench was 
informed, through the space-time definitions, of the maps between the var- 
ious charts. Numerical operations, such as geodesic tracing, are performed 
on a single chart, until a chart boundary or other obstacle is encountered, 
at which point the algorithms are able to transform the data into another 
coordinate system and resume computation there. 

The components of the metric tensor on each coordinate chart, together 
with the maps between charts, define a space-time in GRworkbench. For 
numerical operations which involve derivatives of the metric components, 
such as geodesic tracing, simple, robust numerical methods are employed to 
compute the derivatives. 

A highly general visualisation system was implemented in GRworkbench. 
In the coordinate system of choice, space-times are visualised by transform- 
ing the 4 coordinates under arbitrary distortions down to a 3-dimensional 
visualisation space, which is then rendered on the screen using the OpenGL 
graphics library. Higher-dimensional structures (surfaces, volumes, hyper- 
volumes), such as the event horizon of a black hole, are also intelligently 
visualised under arbitrary distortions. 

FigureEUis a screen-shot from GRworkbench showing a time-like geodesic 
in the Kerr space-time, which describes the gravitational field around a ro- 
tating black hole. The geodesic represents the world-line of a particle falling 
into the near-field of the black hole, orbiting the event horizon several times, 
and then escaping in a different direction. The spherical object in the centre 
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Figure 2.1: GRworkbench screen-shot showing an interesting time-hke 
geodesic in the Kerr rotating black hole space-time. 
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of Figure 12.11 is the event horizon of the black hole. Elements of the GUI are 
visible in the top-left corner. 

The interesting geodesic of Figure 12.11 was obtained in a just a few min- 
utes using the fast turn-around of real-time geodesic tracing and visualisa- 
tion. Other physically interesting situations can be explored visually in a 
similar way. GRworkbench enables users to quickly get an intuitive 'feel' 
for the properties of a space-time, and is thus also potentially useful as an 
educational tool. 

2.3 Objective 

Simple visual experiments have been performed in GRworkbench, demon- 
strating its utility. However, the simulation of more complex physical situa- 
tions was hindered by the numerical methods, which were not as efficient or 
flexible as they could be, and the differential geometric engine, which was not 
sufficiently general for rapid extension. The modification of GRworkbench, 
with the aim of performing complex numerical experiments, is the topic of 
this thesis. 



Chapter 3 

Functional programming 



The numerical and differential geometric engine of GRworkbench has been 
rewritten during 2003 within the framework of functional programming. An 
overview of C++ and functional programming is presented in this chapter. 
The benefits for GRworkbench are discussed in Section [231 Numerical meth- 
ods and differential geometry within this functional framework are the topics 
of Chapters m and [5], respectively. 

3.1 Functions 

In the traditional programming languages of scientific computing, such as C, 
C++, and Fortran, a program typically consists of routines which operate 
on data stored in program variables. Every variable in C++ has a type, 
and there is a natural correspondence between C++ types and standard 
mathematical sets. Table 13.11 lists the most important examples. 

The first two sets in Table 13.11 Z and M, are represented in some way 
or other in every language of scientific computing. The type name double 
stands for 'double-precision floating point number'. 



Set 


C++ type 


Notes 


Z 
M 
M" 


int 
double 

nvector<double> 

function<B (A)> 


max. ±(231 _ 

max. ~ ±10'^°^, precision ~ 10"^^ 
(as for double) 
see Section l3.2l 



Table 3.1: Correspondence between certain sets and C++ types in GRwork- 
bench. 
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The nvector<T> type, written by Antony Searle, uses the C++ template 
mechanism0 to provide a type representing n-tuples of any other type T. The 
type T is called the template parameter. In the case of M", T will be double. 
The template parameter may itself be an n vector, as in nvector<nvector< 
double>>, which is a type representing the set of m x n matrices with real 
entries. 

The following is a routine in C++: 
double mean(double a, double b) 

{ 

return (a + b) / 2; 

} 

The corresponding mathematical definition is 

mean : M x M M, 

mean(a, 6) = — ^— . (3.1) 

The first line of the routine conveys the same information as the first line of 
(13.11) : the routine mean takes two real numbers as arguments, and returns a 
real number. The rest of the routine definition, enclosed in braces, encodes 
the second line of (13. ip . 

The signature of a routine is obtained by taking the first line of a routine 
and removing the routine name and argument names, leaving only their 
types. Thus the signature of the routine mean is double (double, double), 
and the signature of a function /: x Z ^ M would be double (nvector< 
double>, int). 

We may define a function as anything which behaves like the routine mean 
above, in the sense that it accepts zero or more arguments, and returns a 
value. In traditional programming languages (C, Fortran) the only possible 
functions are routines, and so the terms 'function' and 'routine' are used 
interchangeably. The key feature of functional programming is that there 
can be functions other than the routines typed in by the programmer — 
functions created while the program is running. The mechanism to achieve 
this is introduced in Section 13.31 To create functions at run-time, we need 
to be able to store them in variables, which is the topic of the next section. 

3.2 Functions as data 

The capability to store functions in variables is not unique to functional 
programming. Most languages used for scientific computation have some 

iSee [IS], page 327. 
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way to store a reference to a program routine; GRworkbench uses the Boost 
Function Library [6j. The Function Library provides the templatised type 
function<T> representing a function whose signature is T. The following 
code fragment shows how the routine mean can thus be stored in a variableo 

function<double (double, double)> f = mean; 
// the following two lines are now equivalent 
double X = f(l, 2); 
double X = mean(l, 2); 

Observe from the last two lines that the variable f can be used just like the 
routine mean; they are both functions. 

In general, if we let (^i x ■ ■ ■ x An ^ B) denote the set of functions from 
Ai X ■■■ X An to B, then the corresponding C++ type is function<B (Al, . 
. . , An)>, where the sets B, Ai, . . . ,An correspond to the types B, Al, . . . , 
An. The fourth row of Table [3?T] summarises this relationship. 

The most important consequence of the capability to store functions in 
variables is that functions can be arguments to other functions. To illustrate 
this, consider the following routine, which approximates the derivative of a 
function / at a point 

double slope(function<double (double)> f, double x) 

{ 

double h = 0.1; 

return (f(x + h) - f(x - h)) / (2 * h); 

} 

The corresponding mathematical definition is 

slope: (R ^ M) X R ^ M, 

. . f{x + h) - f{x-h) 

slope(/,a;) = — , h = 0.1. (3.2) 

2n 

Again the first line of the routine definition encodes the same information as 
the first line of (13.21) . and the remainder of the routine definition, enclosed 
in braces, encodes the second line of (13. 2p . 

In addition to differentiation, many other numerical algorithms naturally 
take a function as an argument. Two classic examples are 

minimise : (R ^ M) x R ^ R, 

minimise(/, x) = (a local minimum of / near x), (3.3) 

^Anything after the characters // in a Hne of C++ code is a comment, and is ignored 
by the compiler. 

■^This crude method for estimating the derivative is for illustrative purposes only; the 
differentiation algorithm employed by GRworkbench is described in Section [4.31 
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and 

integrate : (M ^ M) x R x R ^ M, 

integrate (/, a, 6) = (numerical estimate of / f{x)dx). (3.4) 

J a 

Finally, note that the signature of slope is double (function<double ( 
double) >, double), and so slope itself may be stored in a variable of type 
function<double (function<double (double)>, double)>. Every function in 
C++ can be stored in a variable of type function<T>, where T is the signa- 
ture of the function. 

3.3 Creating functions at run-time 

Consider the following function, defined in terms of the slope function (13. 2p : 
derivative: (M ^ R) ^ (M ^ R), 

derivative(/) = g, M R, g{x) = slope(/, x). (3.5) 

For any function /, it returns the function which returns the slope of f at its 
argument. 

This expression of the operation of numerical differentiation as a map- 
ping from functions to functions is more flexible than slope. By recursively 
applying derivative, for example, we have derivative(derivative(/)), which is 
an approximation to the second derivative of /. Using only the mechanisms 
introduced so far, however, we cannot encode (13. 5p in C++. 

3.3.1 Functors 

New types are created in C++ by writing a class. A class may optionally 
define an operator(]0 routine, in which case it is called a functor class^ A 
variable whose type is a functor class is a function as defined in Section 13.11 
To see this, consider the following functor class: 

class multiply_functor 

{ 

public: 

// constructor (see the discussion, below) 
multiply_functor(double a_) 

^(pronounced 'operator parenthesis' or 'the parenthesis operator') 

^The use of the term 'functor' in category theory is not related to its use in this context. 
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{ 



a = 



} 



double operator() (double x) 

{ 



return a * x; 



} 



private: 



double a; 



}; 



It can be used in the following way: 

function<double (double)> f = multiply_functor(1.5); 
double y = f(3); 

This code fragment sets f to the function which returns 1.5 times its argu- 
ment, and thus it sets y to 4.5. 

A functor class represents the function encoded by its operator() routine, 
parameterised by the variables in its private: section. The variables in the 
private: section are initialised by the constructor, which always has the same 
name as the functor class. In the code fragment, above, the line a = a_; 
initialises the private: variable a with the value of the variable a_, which was 
passed to the constructor of multiplyJunctor. 

Thus, multiply_functor represents the class of functions which multiply 
their argument by some constant a G M; the value of the parameter a is the 
argument to the constructor. We may even think of the constructor itself as 
a function: 



Using a functor class we can encode the derivative function (13. 5 p in C++: 
class derivative_functor 



multiply Junctor: M ^ (R ^ M), 
multiply_functor(a) = /, /: R ^ M, f{x) = ax. 



(3.6) 




clerivative_functor(function<double (double)> f_) 

{ 



f = f-; 

} 



double operator() (double x) 
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{ 

return slope(f, x); 

} 

private: 

function<double (double)> f; 

}; 

function<double (double)> derivative(function<double (double)> f) 

{ 

return derivative_functor(f); 

} 

If we were to replace the primitive slope routine with a more sophisticated 
algorithm for numerical differentiation, then this derivative routine would be 
a good approximation to the mathematical operation of differentiation. For 
example, derivative(sin) would be a good approximation to the function cos|^ 

3.4 Applicability to GRworkbench 

There are two reasons why functional programming is an ideal framework 
in which to implement the numerical and differential geometric aspects of 
GRworkbench. Functional programming permits numerical operations, like 
derivative, to be expressed in a way which closely resembles the mathematical 
operation that they approximate; and many fundamental notions in differen- 
tial geometry and general relativity, such as the action of the metric tensor, 
and particle world-lines, are functions. 

By elevating functions to the same level as traditional data types (Z, 
M), functional programming makes these notions directly representable as 
variables in C++ code. As we shall see in Chapter [6l this is invaluable in 
the construction of numerical experiments. 



^The functions sin and cos, and many other standard functions, are built-in to C++. 



Chapter 4 
Numerical methods 



The numerical engine of GRworkbench has been rewritten during 2003 within 
the framework of functional programming. Functional algorithms have re- 
placed third-party routines and inline implementations of simpler methods. 
Some algorithms needed to be rewritten or added as part of the develop- 
ment of GRworkbench for numerical experiments, as described in Chapter [HI 
while other changes were directed towards increasing robustness, accuracy, 
or speed of computation. 

A technique for scale-independent computation is described in Section HUTU 
and the method of Richardson extrapolation is introduced in Section 14. 2[ 
These tools are employed in new implementations for the operations of dif- 
ferentiation, integration of ordinary differential equations, and function min- 
imisation, which are described in Sections 14. 3[ 14.41 and 14.51 respectively. 

4.1 Scale-independent computation 

As mentioned in Section [3TTI the name of the type double, which represents 
real numbers in GRworkbench, stands for 'double-precision floating point 
number'. The term 'double-precision' arises from the size of the data type, 
64 bits, being twice that of the smallest floating point data type in C++, 
which is called float and referred to as 'single-precision'. The term 'floating 
point' refers to the particular way that numbers are encoded in the 64 bits. 

Floating point numbers are represented in mantissa- exponent form, which 
is similar to standard scientiflc notation. For example, the number 1.234 x 
10~^^ is represented as a double by 1.234e— 56, where 1.234 is the mantissa, 
which can contain up to 15 signiflcant flgures, and —56 is the exponent, which 
ranges from -308 to +3080 These limitations of the double data type were 

^The mantissa is stored in 52 bits, so its precision is one part in 2^^ ~ 4.5 x 10^^. The 
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summarised in Table I3.1[ 

The alternative to mantissa-exponent form is fixed-point form, in which a 
certain number of bits (32 bits, say) store the part of the number to the left 
of the decimal point, and the remaining bits (31 bits, say) store the part of 
the number to the right of the decimal point, with 1 bit reserved to indicate 
the sign (+ or — ) of the number. In this form, the largest representable 
number is ~ 2^^, and the smallest (in magnitude) representable number is 
~ 2~^^, so the example above, 1.234 x 10~^^, is not representable at all. 
Mantissa-exponent form, offering a wider range of length scales, and the 
same precision at all length scales, is more suitable than fixed-point form for 
scientific computation. 



4.1.1 Approximate equality 

In approximate methods, it is necessary to have a notion of two numbers 
being approximately equal, to some relative precision e. For example, suppose 
e = 0.01; then we want to consider 1.001 x 10*^^ to be approximately equal 
to 1.002 X 10^^, because their difference, 10^°, divided by either of their 
magnitudes, ~ 10^^, is ~ 10"'^ < e. On the other hand, we also want to 
consider to be approximately equal to 10~^, simply because 10"'^ < e. 
We require a definition of approximate equality which satisfies both of these 
examples. 

A notion of approximate equality is also required for elements of other 
sets, most importantly W^, where there is an additional consideration. Con- 
sider two vectors Vi, V2 G M^, 





'10^" 




"10^" 


Vi = 


1 


, V2 = 


2 



(4.1) 



Denoting the standard Euclidean norm on by 
IIV2II 3> 1, ||v2 — vill = 1, and 



we have that ||vi| 



|V2 - vil 

l|vi|| 



< e. 



(4.2) 



However, we may not want to consider the vectors vi and V2 to be approxi- 
mately equal, because their second components are not approximately equal, 
and the scale of interest of the first component may be different to that of 
the second component. 

exponent is stored in 11 bits, so binary exponents up to ±2-'^" — ±1024 can be represented, 
corresponding to decimal exponents of ilogj^g 2^"^^ ~ ±308. 
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In the literature, it is common for numerical algorithms to assume that 
the scale of interest is approximately unity, or at least that it is uniform for 
all components of a vector or matrix; for such algorithms, it is necessary to 
appropriately normalise input variables, and then apply the inverse transfor- 
mation to the output of the algorithm. Definitions like double tiny = l.Oe 
—30; are also common, where the variable tiny is intended to be smaller than 
any quantity that might otherwise arise, apart from zero. Such a definition 
invalidates the routine for scales smaller than 10~^°, which partially nullifies 
one of the main benefits of floating point arithmetic. Whenever either of the 
two issues above was encountered while implementing the numerical meth- 
ods of this chapter, it was found that, by rethinking the relevant parts of the 
algorithm in terms of a general notion of approximate equality, the problem 
could be avoided. 

In the redesigned numerical engine of GRworkbench, the notion of ap- 
proximate equality for any set S is represented by the function 

approx_equal : 5* x 5* x M ^ {true, false}, 

I true, if relative_difference(a, 6) < e; 
approx_equal(a, o, e) = < (4-3) 

I false, otherwise, 

where the function relative_difference encodes, for each set S, a method to 
determine to what precision two given elements are equal. The range of 
approx_equal, {true, false), is represented by the type bool in C++. 

The default definition]! for any set 5* which has a norn|§ || ■ || and is closed 
under an addition operation, is 

relative_difference : S* x S* ^ M, 

relative_difference(a, 6) = — ^ (4.4) 

max(v«M,l) 

Thus, the relative difference is the absolute difference divided by the geomet- 
ric mean of the absolute values, unless the geometric mean is less than unity, 
in which case the relative difference is just the absolute difference. Defini- 
tion (14.41) is not the only conceivable default definition for relative_difference 
that is suitable for M and that is easily generalisable to other sets with 



^The CH — h template mechanism allows for routines which have no particular type 
specified for one or more of their arguments; such a routine may be called with arguments 
of any type for which the routine body makes sense. 

^The norm on M is represented by the function abs, which is built-in to C++. In 
GRworkbench the norm is defined for other types by specialising (overloading) abs to take 
arguments of other types. 
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norms; but it is the definition employed in GRworkbench. The code of the 
relative_difference routine is hsted in Section [A. II 

The relative_difference function is speciahsed for the case S = R", to 
resolve the problem exemplified by (14. 2p : 

relative-difference : M" x M" ^ M, 



relative_difference(a, b) 



\ 



relative_difference(ai, b^y, (4.5) 



where a = (ai, . . . , a„) and b = (6i, . . . , 6„). Thus, the square of the rel- 
ative difference is the sum of the squares of the relative differences of the 
components. 

The specialisation of the relative_difference routine in GRworkbench has 
signature double (nvector<T>, nvector<T>), where T is a template param- 
eter. As such, the componentwise definition (14.51) applies to ra-tuples of any 
set. In particular, recalling that matrices are represented by the type nvector 
<nvector<double>>, by recursively applying (14.51) we find that the square 
of the relative difference of two matrices is just the sum of the squares of the 
relative differences of their components, independent of their representation 
as vectors of vectors. 

More general than the notion of relative difference, as defined in (14.41) and 
(14.51) . is to associate with each set S and norm || ■ || on S" not just a C++ type 
S, representing 5", but also a function of signature double (S), representing 
the norm || ■ ||. The particular norm on S will depend on what the elements 
of S are being used to represent; multiple norms on M", for example, could 
facilitate the correct definition of approximate equality for the two vectors in 
(14.11) , which will depend on the particular meaning of the various components 
of the vectors. 



4.2 Evaluation of limits 

Two of the numerical methods presented in this chapter (that for differenti- 
ation and that for integration of ordinary differential equations) involve an 
algorithm f{h) which approximates the desired solution as a function of a 
small parameter /i G M, such that 

lini/(/i) = (the exact solution), (4.6) 

but such that /(O) is not defined. The limit (14. 6 p must be estimated by 
evaluating f(h) for a finite number of values of h. For very larg^ values of 

^(relative to the scale over which / varies significantly) 
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h, f{h) will be a poor estimate of the limit; but for very small values of h, 
roundoff error in the floating point arithmetic will contribute significantly to 
the value of f{h). 

To see the effect of roundoff error, let 



so that \imh_,Qf{h) = —1 is the derivative of sinx at x = vr. Now, /(O.l) ^ 
—0.998 equals the limit to 2 significant figures, and in general /(10~"), n G N, 
equals the limit to 2n significant figures, if we perform the computation to 
arbitrary precision. However, if we evaluate, say, /(10~^^) using double pre- 
cision floating point numbers, the result is approximately —0.99996, accurate 
to only 4 significant figures. Accuracy is lost because tt + h differs from tt 
only after 12 significant figures jl and so the computed quantity sin(7r + h) is 
only accurate to 4 significant figures. 

4.2.1 Richardson extrapolation 

The purpose of the technique called Richardson extrapolation is to estimate 
the value of the limit (14.61) from several values of f{h), none of which may 
themselves be sufficiently accurate estimates. The basic method is to con- 
struct a polynomial approximation to the function /, and evaluate it at /i = 0. 
That is, evaluate p(0), where p{h) is the unique polynomial of order m fitting 
the m known values {h, f{h)). 

Given the polynomial of order m passing through m known values, it is 
possible to efficiently determine the polynomial of order m+1 passing through 
the m+1 points consisting of the m original points plus one additional point. 
As such, if the estimate of the limit fl4.6l) afforded by the first m evaluations 
of f{h) is not sufficiently accurate, another single evaluation can be made 
and a new estimate of the limit obtained. 

If the estimate after m+1 function evaluations is approximately equal 
to the estimate after m function evaluations, to within the desired relative 
precision e, in the sense defined in Section 14.1.11 then no more function 
evaluations are made. The most recent estimate, namely the estimate after 
m+1 function evaluations, is then the output of the Richardson extrapolation 
process: an approximation of the limit (14. 6p . 

Richardson extrapolation is particularly useful when a power series of the 
function f{h) about h = is known to contain only even powers of h; this 
is the case for both of the applications of Richardson extrapolation in this 



sin {tt + h) — sin tt sin (tt + h) 



(4.7) 



h h 



^(out of the 15 or at most 16 significant figures representable in the double data type) 
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chapter. The power series may then be treated as a polynomial in h'^, rather 
than a polynomial in h. The extrapolation polynomial is then p{h'^), passing 
through known values {h'^,f{h)). In evaluating the function / at, say, half 
the previous value of h, a new polynomial fitting point is obtained which is 
four times closer to zero. 

In GRworkbench, the templatised class richardson_extrapolation<T>, whose 
code is listed in Section [Al2l represents the operation of Richardson extrapo- 
lation on a function from M to the set represented by the type T; typically T 
is double or an nvector type. The refine routine of the richardson_extrapolation 
class takes one argument of type double and one argument of type T, repre- 
senting a new known value pair {h, f{h)); using the new values, and the values 
supplied in previous calls to the routine, refine computes a new estimate of 
the limit (14. 6p . and computes the difference between the new estimate and 
the previous estimate as an approximation of the error. The most recent 
estimate and error are accessed, respectively, through the routines limit and 
error of the richardson_extrapolation class. 



4.3 Differentiation 

Numerical differentiation in GRworkbench is implemented in terms of the 
class richardson_extrapolation of Section |4. 2. 1| exposing a functional interface 
similar to that developed for the derivative routine of Section For a vector 
space V, numerical differentiation is encoded in a routine 

derivative: (M ^ V) x M x R ^ (M -> V), 
derivative(/, ^,e) = g, g: M. ^ V, 

g{x) = (the derivative of / at x, to relative precision e), (4.8) 

where the argument /i is a characteristic length scale over which the function 
/ varies significantly. Depending on the choice of fi, the routine may not 
successfully converge to an estimate of the derivative to relative precision e. 
The code of the derivative routine is listed in Section IA.3I 

The function g in (14.81) employs Richardson extrapolation to estimate the 
value of 

lin.l^^±^^i^-li^ = limdih), (4.9) 

which is the centred difference approximation to the derivative of / at x. 
Observe that d{h) is an even function of h\ hence a power series expansion 
of d{h) about h = contains only even powers of h, and the Richardson ex- 
trapolation can be performed using the value pairs {h'^, d{h)), rather than the 
value pairs {h, d{h)), with the advantage described at the end of Section l4.2.1[ 



4.3. DIFFERENTIATION 



21 



The first value of h for whicli d{h) is computed by the derivative routine 
is h = fi, the characteristic length scale of the function /; the nth value 
of h is /i/cr"~^, where a = 1.7 is a constant parameter of the algorithm. 
At most nmax = 13 values of h are processed, after which the algorithm 
terminates, and the derivative is undefined. Thus, the algorithm explores 
the region around x at length scales between fi/a""^''^ ~ 10~^/i and fi. The 
particular values of the constants a and rimax were empirically chosen to 
optimise computation speed for the applications of GRworkbench discussed 
in this thesis. 

Previously in GRworkbench, numerical differentiation was accomplished 
by, where an algorithm required it, evaluating d{h) at progressively smaller 
values of h, until the difference between two successive evaluations was smaller 
than the desired precision. The new implementation, employing Richardson 
extrapolation and the C++ template mechanism, converges faster and more 
accurately, and its interface is more general, in that functions from M to any 
sensible set can be differentiated. 

4.3.1 Gradient 

The gradient of a function of M" is defined in terms of derivative. For any 
vector space V, the gradient is defined by 

gradient: {W V) ^ (M" ^ 1/"), 
gradient^) = (7, g-.R^'^V^ 

g{x.) = (derivatives of / at x with respect to the n components). (4.10) 

The code of the gradient routine is listed in Section PV.3.11 

Like many routines in GRworkbench that employ derivative, gradient uses 
default values of /i = 1 and e = 10~^ for the arguments to derivative. In 
general, these routines should be extended to accept these parameters as 
arguments, and to pass them on to all numerical routines which require 
them; the scale information /i in GRworkbench must originally be supplied 
with definitions of the metric. For current applications, the metrics input to 
GRworkbench have unity as an appropriate length scale, and so this extension 
has not yet been performed. 

Previously in GRworkbench, the gradient of a field was computed by, 
where an algorithm required it, explicitly calculating the numerical deriva- 
tives with respect to the various components of the vector argument, and 
populating a vector with the results. Like derivative, the new implementation 
employs the C++ template mechanism to create a more general algorithm, 
which can apply the definition f l4.10p for any set V for which it makes sense. 
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4.4 Integration of ordinary differential equa- 
tions 

Previously in GRworkbench, numerical integration of ordinary differential 
equations (odes) was performed using the third-party Slatec ddrivS Runge- 
Kutta algorithm [16j, originally written in Fortran, converted to C using a 
Fortran-to-C source code converter, and then adapted to the C++ code of 
GRworkbench. During the course of the project, it was discovered that the 
Slatec algorithm was coded such that only one numerical integration can 
be in operation at any time. Normally, this presents no problem; but in 
the case that the function / which gives the derivatives in the initial value 
problem specification. 



is defined in terms of the integration of another, separate ODE, the Slatec 
algorithm is inadequate. 

It was decided that, rather than further adapting the Slatec algorithm, 
a general ODE integrator should be directly implemented in the newly func- 
tional framework of GRworkbench. The Bulirsch-Stoer method, described in 
[T2] . pages 724-732, and pages 484-486, was selected based on argu- 
ments in [18], pages 487-488, which recommend it for odes whose derivative 
functions / are smooth]^ and for applications where high accuracy is required. 
The Bulirsch-Stoer method is generally inferior to Runge-Kutta methods for 
ODES for which the derivative function / contains discontinuities near the ex- 
act solutionjzl or for stiff ODEs, but neither of these cases occur in the current 
applications of GRworkbench. 

The Bulirsch-Stoer method, as implemented in GRworkbench., applies 
Richardson extrapolation to a series of estimates obtained using the modified 
midpoint method, from [12], pages 722-724. The modified midpoint method 
is an algorithm for estimating y{H) from y(0), evaluating the derivative 
function / at the initial point yo and at n other points, by the following 



^By 'smooth' we mean not varying significantly on scales much smaller than the region 
of integration. 

^(because Bulirsch-Stoer steps are longer than Rungc-Kutta steps, and are thus more 
likely to 'accidentally' land on or near a discontinuity) 




y(0) = yo, 



(4.11) 



4.4. INTEGRATION OF ORDINARY DIFFERENTIAL EQUATIONS 23 



process: 

h = H/n, 
zo = yo, 

Zi = Zo + /l/(0,Zo), 

Zm+i = Zm_i + 2hf{mh, z.^), 

y(if)^^!i±^!^zi±M^lM. (4.12) 
2 

It is a second-order method in h. 

The modified midpoint estimate of y{H), for the initial value problem 
(14.111) . is a function m{h). The modified midpoint method is chosen for 
extrapolation using Bulirsch-Stoer because, like the function d{h) in (14. 9 p 
employed by derivative, in a power series of m{h) about h = 0, all odd powers 
of h cancel out, and so the extrapolation can be performed in h"^. 

The modified midpoint method is represented in GRworkbench by the 
class modified_midpoint_stepper, whose code is hsted in Section [A. 4. 1[ It must 
be supplied with the derivatives function / and the initial data yo. The only 
routine, step, takes H and n as arguments, and returns the estimate y{H). 

The difficult problem of choosing the optimal value for H, so that the 
Richardson extrapolation will not take too many steps, but so that a signif- 
icant distance in x will be covered, is discussed in [12], pages 726-728. 

The class bulirsch_stoer, whose code is hsted in Section lA.41 is adapted 
from the implementation of the Bulirsch-Stoer method in [12]. The class 
must be supplied with the same information as modified_midpoint_stepper, 
as well as: a characteristic length scale in x, over which / in (14.111) varies 
significantly; the maximum number of step^to try before giving up; and the 
desired relative accuracy of the solution. The routine step takes an argument 
indicating the desired final value of x, after which the routines x and y return, 
respectively, the final values of x and y obtained by the algorithm; if the 
result of a call to the routine x equals the argument given to step, then the 
integration was successful. 

The new implementation of numerical ODE integration in GRworkbench 
is more general than the Slatec Runge-Kutta algorithm. Previously, the 
ODE integrator required the function / to satisfy /: x R ^ R", and to be 
encoded using the built-in array notation of C++ (rather than in terms of 
nvector, or some other type). Now, the function / can satisfy / : V xM. V , 
where V is any vector space. 



*A step is a successful Richardson extrapolation of the results of as many calls to 
modified_midpoint_stepper as are necessary. 
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4.5 Minimisation of functions 

Previously, the applications of GRworkbench did not necessitate a mechanism 
to find local minima of functions. The development of tools for numerical 
experimentation, as described in Chapter [6l highlighted the need for a general 
algorithm which, for a function / : R" —>■ M, can locate a minimum of / near 
a given initial 'guess' point x. 

If / : R ^ M, then a local minimum of / can be bracketed by three num- 
bers a < b < c which satisfy /(a) > f{b) < f{c). More efficient algorithms 
exist for this special case; GRworkbench employs Brent's method, from [T2], 
pages 402-405, which repeatedly refines the bracket on a minimum by fitting 
the three smallest function values found so far (the smallest of which will 
be /(&)) to a parabola, and using the exact minimum of that parabola as 
the next trial point; it converges quadratically near the minimum. Brent's 
method is represented in GRworkbench by the functor class brent_minimiser, 
whose constructor must be supplied with the function /; it is then the func- 
tion 

brent _minimiser :MxMxM^]RxM, 
brent_minimiser(xo, /U, e) = (x^in, fix^nm)), (4.13) 

where x^i^ is within relative precision e of a local minimum of / near xq, and 
/i is a characteristic length scale over which / varies significantly. The code 
of the brent_minimiser class is listed in Section [A. 5. 11 

4.5.1 Multi-dimensional minimisation 

In the general case of multi-dimensional minimisation, minima cannot be 
bracketed, and minimisation consists, more or less, of 'rolling' downhill from 
the initial guess xq. GRworkbench employs Powell's method, from [12], pages 
412-418, which proceeds by using brent_minimiser to minimise the function 
one-dimensionally in each of n linearly independent directions. The n basis 
directions are then updated, based on the overall distance moved from xq, 
and the process is repeated with the new directions. The problem of how to 
choose the right basis directions is discussed in [12]. 

Powell's method is represented in GRworkbench by the functor class 
powelLminimiser, whose constructor must be supplied with the function / : M" 
M; it is then the function 

powelLminimiser: M" x M„xn x M ^ M" x M, 
powelLminimiser (xo, B, e) = (x^in, /(xmin)), (4.14) 
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where ^ram is within relative precision e (in the Euchdean norm on M") of a 
local minimum of / near xq, M„xn is the set of n x n matrices with real entries, 
and B is the matrix whose columns are the initial directions to minimise over. 
The minimisation is made over the subspace of R" spanned by the columns 
of B] this will be all of R" only if the columns of B are linearly independent. 

The code of the powelLminimiser class is listed in Section |A.5[ The imple- 
mentation of Powell's method in [12] requires a separately coded implemen- 
tation of Brent's methocl^ to perform the minimisations over one- dimensional 
subspaces of M"; the quite general interface of the brent_minimiser class makes 
this inelegance unnecessary in the implementation of Powell's method in GR- 
workbench. 

4.6 Conclusion 

The rewritten and extended numerical engine of GRworkbench is more effi- 
cient, robust, and general. The implementation of sophisticated algorithms 
for key operations yields increased computation speed. The relative_difference 
abstraction enables algorithms to be encoded with consistent notions of ap- 
proximate equality, making them more robust and elegant. Through the 
C++ template mechanism, numerical methods can be encoded such that 
they can be applied to any sets which have the required structure defined 
upon them. 



^See [E], pages 418-419. 



Chapter 5 

Functional differential 
geometry 

The differential geometric engine of GRworkbench lias been rewritten within 
the framework of functional programming, using the functional numerical 
tools of Chapter HI The definition of charts, and the components of the 
metric on charts, is discussed in Section 15.11 Collections of charts, and inter- 
chart maps, are introduced in Section 15.21 The representation of points and 
tangent vectors as C++ classes is described in Section [531 

Table 15.11 summarises the correspondence between important concepts 
in differential geometry and their representations in GRworkbench. Each 
correspondence is described in detail in this chapter, but, as the concepts are 
interrelated. Table 15.11 will be useful when reading the earlier sections. 

5.1 Charts and the metric components 

A chart is a subset C C M", representing a coordinate system on a subset 
Aic C of the space-time manifold Ai. We denote by (pc- -Mc ^ C the 
one-to-one and onto function which maps points in A^c* into the chart C. 

A space-time in GRworkbench consists of the definition of the components 
of the metric tensor on one or more charts, and the definition of maps (coor- 
dinate transformations) between those charts. In this section we describe the 
definition of the metric components on charts; discussion of the inter-chart 
maps is deferred until Section 15.21 

The coordinates of a point on a chart, {x*}"^]^ G M", or simply G M", 
where n is the dimensionality of the space-time, are represented by a variable 
of type nvector<double> (see Table 13. ip . The components of the metric 
tensor Qab at a point on a chart are represented as an n x n matrix, by a 
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Concept 


Representation in GRworkbench 


Section 


Space-time 


atlas 




Coordinates 


nvector<double> 




Metric components 


nvector<nvector<double>> 




Inter-chart map 


See fl5.9l) 


15.21 


Point 


point 


15.4.11 


Tangent vector 


tangent_vector 


15.4.31 


Metric 


function<double (tangent_vector, tangent_vector)> 


15.4.41 


World-line 


function<point (double)> 


15.4.21 



Table 5.1: Representation of important differential geometric concepts in 
GRworkbench 



variable of type nvector<nvector<double>>. A function which defines the 
metric components Qab, as a function of the chart coordinates x*, might then 
be of the form 

chart: M" M^xn, 

chart(x*) = gab\x^, (5.1) 

represented in GRworkbench by a function of signature nvector<nvector< 
double>> (nvector<double>). In general, however, the chart coordinates 
are an open subset of M", and so (15.11) will not be defined everywhere in M". 
A mechanism is required to represent functions which are only defined on a 
subset of some other, standard, setQ 

5.1.1 The optional mechanism 

GRworkbench employs the Boost Optional Library ^ to represent functions 
which are undefined for some values of their arguments. The Optional Li- 
brary provides a templatised type optional<T>, which represents the set 
S U {0}, where S is the set corresponding to the template parameter type 
T, and is a special value taken by functions at points where they are 
undefined. 

The optional template might be used in the following way: 

optional<double> square_root(double x) 

{ 

if (x < 0) 

^By 'standard set' we mean a set which is represented by a type in C++, such as those 
in Tables O and O 
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{ 



// undefined; return the special value 'undefined' 
return optional<double>(); 



} 



else 



{ 



// defined; return the result of the standard C++ square root 

algorithm, 'sqrt' 
return optional<double>(sqrt(x)); 



} 



} 



Thus, by returning a variable of type optional<double>, instead of a variable 
of type double, the square_root routine can return the special value (using 
the code return optional<double>();) to indicate points where the algorithm 
is undefined; in this case, is returned for negative values of the argument 



The optional mechanism is most useful when the caller of a function can- 
not know beforehand whether the function will be defined at the arguments 
to be given to it. This would be the case for callers of the function (15.11) : 
the differential geometric algorithms in GRworkbench must be coded in such 
a way that they can operate on any space-time definition, without prior 
knowledge of the particular coordinate systems (charts) they will be working 
in. 

We can now modify (15.11) to support charts defined on subsets of M", using 
the optional mechanism. Thus, in GRworkbench, functions which return the 
metric components Qat, as a function of the chart coordinates x^, are of the 
form 



The corresponding C++ type is 

function<optional<nvector<nvector<double>>> (nvector<double>)>, (5.3) 

for which GRworkbench declares a short synonym, chart, using the C++ 
typedef mechanism: 

typedef function<optional<nvector<nvector<double>>> (nvector<double 

>)> chart; 

References to charts are stored in variables of type shared_ptr<chart>, using 
the Boost Smart Pointers Library 



X. 



chart: M' 




n 



M„xn U {0}, 

gab\x', if the are valid chart coordinates; 
0, otherwise. 



(5.2) 
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5.1.2 Example chart and metric components 

In this section we demonstrate the encoding of the flat space metric of special 
relativity, in cylindrical coordinates {t, r, (j), z), using a C++ function of type 
chart. The line element is 

ds'^ = -dt^ + dr^ + d^'^ + dz'^ , (5.4) 

so the metric components, as functions of the chart coordinates {t, r, 0, z), are 
~gtt = Qrr = dzz = ^, g<t><t> = r"^, and all other gab = 0. The chart coordinates 
are valid in the open subset of M" satisfying 

t G ( — CXD, oo), 
r e (0,oo), 
(0,27r), 

z G (— oo, oo). (5.5) 

The following routine encodes (15.41) and (15. 5p in C++: 

optional<nvector<nvector<double>>> flat_metric_cylindrical(nvector< 
double> x) 

{ 

// t = x[0], r = x[l], phi = x[2], and z = x[3] 
if (x[l] <= or x[2] <= or x[2] >= 2 * pi) 
{ 

// invalid chart coordinates; return 'undefined' 
return optional<nvector<nvector<double>>>(); 

} 

else 

{ 

// valid chart coordinates; compute and return metric components 
nvector<nvector<double>> gab; 

gab[0][0] = -1; 
gab[l][l] = 1; 
gab[2][2] =x[l] *x[l]; 
gab[3][3] = 1; 
// all other gab = 



} 

} 



return optiona I < nvector< n vector<double> > > (ga b) ; 



The operator [i], applied to an nvector such as in x[i], returns the ith compo- 
nent of the vector. 
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The opening if statement determines whether the argument x represents 
vahd chart coordinates; if so, the metric components are computed in the 
variable gab, and returned; if not, is returned. All space-times in GR- 
workbench have the metric components defined on each of their charts by 
functions like flat_metric_cylindrical, above. 

5.1.3 The connection 

The components of the connection, or the Christoffel symbols, are the useful 
quantities defined in terms of the metric components Qab by 



where gah,c denotes partial differentiation of Qah with respect to the coordinate 
x'^, and g""^ denotes the contravariant components of the metric tensor. The 
Christoffel symbols are used by the numerical differential geometric functions 
of Chapter El 

The GRworkbench routine connection accepts an argument of type chart, 
and returns a variable of type function<optional<nvector<nvector<nvector< 
double>>>> (nvector<double>)>, representing the function which returns 
the components (15. 6p as a function of the chart coordinates. 

The differentiation of the metric components Qab is accomplished using 
the numerical tools of Chapter HI A function which returns the components 
of gab,c, as a function of the chart coordinates, is given simply by gradient(c), 
where c is the function, of type chart, which returns the metric components 
gab as a function of the chart coordinates. 

The matrix of contravariant components g"''^ of the the metric is simply 
the inverse of the matrix g^b of covariant components. This matrix inversion 
is performed in GRworkbench using standard row reduction techniques (see 
for example 0, pages 115-116). 

5.2 Inter-chart maps 

As mentioned at the beginning of Section 15. H space-times are defined by 
specifying, together with the metric components on each chart, maps between 
the various charts. 

For two charts A,Bg M", the inter-chart map from A to i? is 



^ab — 7,9 '^{gad,b + gbd,a — gab,d), 



(5.6) 



(f)AB- A^ B, 



(5.7) 
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where 4>b\ma function 0^ restricted to the set M.a, and o denotes 

function composition. The inter-chart maps must be specified to complete 
the definition of a space-time. 

In the definition (15. 7p . the domain A of (j)AB is, in general, a subset 
of M". Hence (pAB cannot be represented by a variable of type function< 
nvector<double> (nvector<double>)>; instead, the optional mechanism of 
Section 15.1.11 is again employed. Thus, in GRworkbench, an inter-chart map 
from a chart A to a chart B is represented by a function of the form 



map: ^ M" U {0}, 



/ »x _ i JbIma ° if (x') e A and 0^ (x*) G Mb; ,^ 

mapi^x J — ^ • '^'^■"^ 

0, otherwise. 

The corresponding C++ type is 

function<optional<nvector<double>> (nvector<double>)>. (5.9) 

As with charts, the C++ typedef mechanism is used to define a synonym 
map for the type (15.91) . References to maps are stored in variables of type 
shared_ptr<map>. 



5.2.1 Example inter-chart map 

In this section we demonstrate the encoding in GRworkbench of an inter-chart 
map of the form (15.91) . which transforms between two cylindrical coordinate 
systems like example (15.51) in Section 15.1.21 with the coordinate systems dis- 
placed from each other by vr in the (f) coordinate. Together, the two coordi- 
nate systems thus cover the entire flat-space manifold of special relativity, 
except for the line r = 0. 

The coordinate transformation, of the form (15.81) . is 

revolve: M" ^ U {0}, 

{(t,r,0 + 7r,2;), if < tt; 
0, if = vr; (5.10) 

(t, r, — TT, 2;), otherwise, 

and is encoded in C++ in the following way: 

optional<nvector<double>> revolve(nvector<double> x) 

{ 

// t = x[0], r = x[l], phi = x[2], and z = x[3] 
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if (x[2] == pi) 

{ 

// mapping not defined; return 'undefined' 
return optional<nvector<double>>(); 

} 

else 

{ 

// mapping defined; perform transformation 
nvector<double> y; 

y[0] = x[0]; 

y[i] = x[l]; 
if (x[2] < pi) 

y[2] = x[2] + pi; 

else 

y[2] = x[2] - pi; 
y[3]=x[3]; 

return optional<nvector<double>>(y); 

} 

} 

The operator ==, used in the first if statement, is the test for equality in 
C++. 

By using a functor class (Section 13.3. II) . we could parameterise the trans- 
formation revolve on the angle of rotation, which is currently vr. All space- 
times in GRworkbench have their inter-chart maps specified by routines or 
functors like revolve, above. 

5.3 Atlases 

A collection of charts with the metric components defined on them, of the 
form (15. 2p . and a collection of inter-chart maps, of the form (15. 8p . together 
comprising a space-time, are represented in GRworkbench by the class atlas. 
The atlas class uses C++ Standard Template Library (stl) ^lOj containers 
to maintain the collections of charts and maps. 

An atlas contains a std::set of charts, and a std::map from std::pairs of 
charts to inter-chart map definitions of type mapl^ An atlas also contains an 
int named dimension which stores the dimensionality of the space-time. 

The members charts and maps of class atlas are used by the differential 
geometric algorithms of GRworkbench to, respectively, enumerate the set of 



^std::set, std::map, and std::pair are STL templates; see [lOj . 
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all charts, and retrieve the inter-chart map between any two charts. If two 
charts do not overlap at all, there will be no inter-chart map between them; 
this is equivalent to there being an inter-chart map between them that always 
returns 0. 

5.4 Points and tangent vectors 

For a point, a valid chart is a chart containing the point; for a tangent vector, 
a valid chart is a chart containing the point whose tangent space contains 
the tangent vector. While points and tangent vectors may be represented 
by their coordinates on a valid chart, it is useful to have a representation of 
these objects which is not linked to any particular chart. The GRworkbench 
representation for points is described in Section [5.4.H and the representation 
for tangent vectors is described in Section 15.4.31 

5.4.1 Points 

The abstract notion of a point p & Ai, independent of any particular coor- 
dinate system, is represented in GRworkbench by the class point. A point is 
constructed from three pieces of information: the atlas to which it belongs, 
a chart which contains it, and its coordinates on that chart. 

The context and valid_chart routines of class point return, respectively, the 
atlas and the chart from which the point was constructed. Numerical opera- 
tions involving points can only be performed in terms of a valid coordinate 
system, so the valid_chart routine is used whenever a variable of type point is 
an argument to a numerical differential geometric routine in GRworkbench. 

Change of coordinates 

The operator[] routine of class point, which takes one argument, a variable of 
type chart, returns a variable of type optional<nvector<double>>, represent- 
ing the coordinates of the point on the given chart. (The optional mechanism 
of Section 15.1.11 is used because a particular point may, or may not, have 
coordinates on the given chart.) Thus, if p is a variable of type point, and c 
is a variable of type chart, then the coordinates of p on c are given by p[c]. 

Let a be the variable of type chart from which p was constructed. If c and 
a represent the same chart, then p[c] will simply return the coordinates from 
which p was constructed. If, on the other hand, c and a are different charts, 
then GRworkbench will use the maps member of the atlas class to determine 
if there is an inter-chart map from a to c defined; if so, then the inter-chart 
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map is used to compute the coordinates of p on c, which are then returned; 
if not, then is returned, indicating that p is not contained in the chart c. 

5.4.2 World-lines 

A curve in space-time, such as a world-hne, is a function A: M Ai; 
such functions are represented by variables of type function<point (double) 
>. However, if the curve A is not defined for all values of its real parameter, 
then it will instead be represented by a variable of type function<optional 
<point> (double)>. All curves in GRworkbench are in fact represented in 
this latter form, because they are often defined in terms of numerical pro- 
cesses which may not converge to a solution. The computation of geodesies, 
discussed in Section 16. 2[ exemplifies this. 

The C++ typedef mechanism is used to define the synonym worldline for 
the type function<optional<point> (double)>: 

I typedef function<optional<point> (double)> worldline; 

5.4.3 Tangent vectors 

The abstract notion of a tangent vector v e Tp, where Tp is the tangent space 
of a point p & A4, is represented in GRworkbench by the class tangent_vector. 
Like a point, a tangent_vector is constructed from three pieces of information: 
the point to whose tangent space it belongs, a chart containing that point, 
and the contravariant component^ of the tangent vector on that chart. 

The context routine of class tangent_vector returns the point from which 
the tangent vector was constructed; through the valid_chart routine of this 
point, a valid chart for the tangent vector can be obtained. As with the point 
class, the operator[] routine of the tangent_vector class, taking one argument, 
a variable of type chart, returns the components of the tangent vector on the 
given chart, in a variable of type optional<nvector<double>>. 

Change of coordinates 

As with the point class, when the components of a tangent vector are re- 
quested on a chart other than that from which the tangent vector was con- 
structed, GRworkbench uses the inter-chart map, if it exists, to compute the 
components. If f * are the components of a tangent vector f at a point p 



■^Whenever we discuss the components of a tangent vector, we always mean its con- 
travariant components. 
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on a chart with coordinates x\ then the components on another chart, with 
coordinates , are 



v' = A\v\ (5.11) 

p 

The columns of the matrix are the derivatives of the inter-chart map 
(p: M" M" with respect to the coordinates of its argument, evaluated 
at p. GRworkbench computes A] , and thereby the components f * , by using 
the methods of Chapter H] to numerically evaluate the derivatives. 

5.4.4 Tangent vectors and the metric 

At a point p, the metric gab is naturally considered as the inner product 

metric : Tp x Tp — R, 

metric(M, f) = Qabu'^v'^. (5.12) 

If M = f in (15.121) . then the sign of metric(M,-u) determines whether u is 
space-like, null, or time-like. If metric(M,M) = —1 then u represents the time 
direction of a physical observer — this is discussed in Section I6.1.1[ 

The function (15.121) is encoded in GRworkbench in the routine metric, 
whose signature is double (tangent_vector, tangent_vector). Also, the operator 
* routine of the class tangent_vector is defined to call metric, so that if u and 
V are variables of type tangent_vector, then the expression u * v is equiva- 
lent to the expression metric(u, v). This notation is reminiscent of the two 
equivalent forms 

gabu'^v'' = Ubv' (5.13) 
for the inner product of two vectors. 



5.5 Conclusion 

The implementation of the differential geometric structure of GRworkbench 
within the framework of functional programming, using the numerical meth- 
ods of Chapter m is robust and elegant. The representation of abstract ob- 
jects such as points and tangent vectors, independent of any particular chart, 
will be useful in the construction of the numerical experiments of Chapter [6l 



Chapter 6 

Numerical experiments 



A numerical experiment is a model of a physical situation in GRworkbench, 
from which a measurement of a physical quantity is obtained. Tools for simu- 
lating physical situations in GRworkbench have been implemented using the 
methods of Chapters H] and [51 Basic operations on points and tangent vec- 
tors are described in Section [UTTl Geodesic tracing and the parallel transport 
operation are the topics of Sections 16.21 and 16.31 respectively. In Section 16.41 
we discuss methods for finding geodesies that are defined implicitly in terms 
of boundary conditions. 

In Chapter [HI the methods of this chapter are used to numerically inves- 
tigate the claim to be discussed in Chapter [71 



6.1 Basic operations 

In this section we describe some operations on points, tangent vectors, and 
world-lines, which will be useful for constructing numerical experiments. 



6.1.1 Tangent vectors and observers 

As was mentioned at the end of Section 15.4.31 a tangent vector u, such that 
metric(M, u) = —1, represents the proper time direction of a physical observer. 
More precisely: physical observers are defined by their time-like world-lines, 
with parameter t; if the tangent vector u to the world-line always satisfies 
metric('U,M) = —1, then the parameter t is the (proper) time coordinate in 
the frame of reference of the observer. 
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Normalisation of a tangent vector is defined by 



normalise : Tp ^ T^ 
normalise (m) = — ^ 



^y\metnc{u, u)\ 



u 



(6.1) 



Thus, the normalisation of a vector m is a vector v such that metric(f,f) = 
±1, according as whether u was space-like or time-like. The definition fl6.1l) 
is encoded in GRworkbench in the routine normalise, which has signature 
tangent_vector (tangent_vector). 

Also useful is the operation of orthonormalisation. The orthonormalisa- 
tion of a vector u with respect to another vector v is defined by 

orthonormalise : Tp x Tp ^ Tp, 

orthonormalise(ti, v) = normalise(metric(u, v)v — metric(f , v)u), (6.2) 

which is encoded in GRworkbench in the routine orthonormalise, which has 
signature tangent_vector (tangent_vector, tangent_vector). Orthonormahsation 
has the property that, if w = orthonormalise(ti, f ), then metric(f,w) = 0, 
and either metnc{w,w) = 1 or metnc{w,w) = —1. 

6.1.2 Orthonormal tangent bases 

An orthonormal tangent basis for Tp at a point p is a set of n vectors in 
Tp that are mutually orthonormal. The metric components expressed in 
an orthonormal tangent basis form a diagonal matrix; this will be useful in 
Section The determination of an orthonormal tangent basis is also called 
diagonalising the metric. 

GRworkbench constructs an orthonormal tangent basis by finding the 
eigenvectors of the matrix g of metric components gab- The eigenvectors are 
orthogonal, because the matrix g is symmetric. The process of determining 
the eigenvectors of a matrix is represented in GRworkbench by the class 
eigen, which is constructed from a variable of type nvector<nvector<double 
>>, representing the matrix whose eigenvectors are to be determined. The 
routine vectors of class eigen then returns a variable of type nvector<n vector 
<double>>, representing the n eigenvectors, and the routine values of class 
eigen returns a variable of type nvector<double>, a list of the corresponding 
eigenvalues. 

The eigen class uses an iterative method to find the eigenvectors of g (see 
p!5] . page 25). Starting with a coordinate basis vector ei, the sequence of 
vectors g'^ei converges, as n — cxd, to an eigenvector vi of g. A second eigen- 
vector t>2 is obtained by seeding the process with 62. Because the sequence 
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g^Ci will tend to converge to the eigenvector which has the largest eigenvalue, 
each successive estimate is orthogonalised with respect to the previously de- 
termined eigenvectors, before the next left-multiplication by g. Once this 
process has been completed, starting with each coordinate basis vector, the 
full set of orthogonal eigenvectors are known. 

If the metric is Lorentzian, then one of the eigenvectors will have a nega- 
tive eigenvalue, corresponding to a time-like direction, and all the others will 
have positive eigenvalues, corresponding to space-like directions. The nor- 
malised eigenvectors constitute an orthonormal tangent basis. The GRwork- 
bench routine orthonormaLtangent_basis takes one argument of type point, 
and one argument of type chart, and uses the eigen class to return a variable 
of type nvector<nvector<double>>, representing a matrix whose columns 
are the components of an orthonormal basis of the tangent space of the given 
point in the given chart. 

6.1.3 Coordinate lines 

If a particular coordinate system on a space-time has known properties, such 
as the metric being independent of one of the coordinates, then it may be 
useful to specify space-time curves explicitly in terms of the coordinates. 
Straight lines in a particular coordinate system are obtained in GRworkbench 
through the coordinate_line routine, which takes three arguments: a point on 
the curve; the chart on which the curve is to be a straight line; and an nvector 
<double> giving the components of the tangent vector to the coordinate line 
at the given point. 

The coordinateJine routine returns a variable of type worldline, as defined 
in Section I5.4.2[ If the coordinate line intersects a chart boundary, then it is 
undefined beyond it; hence the use of the optional mechanism. 

6.2 Geodesies 

Geodesies, the straightest possible lines in a curved space-time, are physically 
important. Geodesies whose tangent vectors are time-like are the world- 
lines of freely-falling observers; geodesies whose tangent vectors are space- 
like represent straight 'rulers', for observers whose world-lines intersect them 
orthogonally; and geodesies whose tangent vectors are null represent the 
world-lines of photons. 

Geodesies are uniquely defined by a point p on the geodesic and the 
tangent vector v of the geodesic at p. The coordinates x'^ of a geodesic on a 
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chart A, as functions of an ajfine parameter t, satisfy the geodesic equation, 

cPx^ „^ dx"" dx'^ 



which involves the connection (15.61) . Note that the components of F^^ in (16.31) 
are a function of the coordinates x'^. 

The equation fl6.3l) is a system of n second order odes in the coordinates 
x^; we may rewrite it as a system of 2n first order odes. Together with 
the n components of an initial point p on A, and the n components of an 
initial vector f G Tp on A, (16.31) defines an initial value problem, which can 
be solved on the chart A using the numerical ODE integration techniques of 
Section lOl 

In general, no single chart will cover the entire space-time. Equation 
(16.31) can only be integrated up to a chart boundary; beyond that, the metric 
components Qab, and hence the Christoffel symbols V^j^, are undefined on that 
chart. 

Let y be a point near the boundary of a chart A, beyond which numerical 
integration of (16.31) fails. If there is another chart B containing y, and an 
inter-chart map from A to B, then integration of (16.31) can be attempted on 
B: Using the inter-chart map, the components x^, in (16. 3p . can be computed 
on B from those on A; using (15. lip , the components dx^/dt, in (16. 3p . can be 
computed on B from those on A; and, using (15. 6p . the components of F^^ at 
y can be computed on B. 



6.2.1 Implementation in GRworkbench 

A point on a geodesic, and the tangent vector to the geodesic at that point, 
are represented in GRworkbench by a variable of type tangent_vector. (The 
context routine of class tangent_vector returns the point at which the tangent 
vector exists.) To determine a new tangent_vector on the geodesic, at a desired 
value t = tfinai of the affine parameter, GRworkbench uses the operator[] 
routines of the classes point and tangent_vector to obtain the initial data for 
equation (16. 3p on each chart, one by one, until it finds a chart on which (16.31) 
can be integrated. 

If no chart exists on which (16.30 could be successfully integrated to the 
desired value tgnai of the affine parameter, then integration to affine parameter 
ifinai/2 is attempted, followed by integration to affine parameter tfinai- If 
either of these integrations fail, then the corresponding interval in t is further 
subdivided, up to a maximum of 7 bisectionso If the maximum number of 

"'^The maximum number of bisections, 7, was empirically determined to be adequate for 
current applications of GRworkbench. 
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bisections is reached without successful integration to t = tfinai, then is 
returned, indicating that the geodesic is undefined at the value tfinai of the 
affine parameter. 

This definition of a geodesic from its initial data is represented in GR- 
workbench by the functor class geodesic, which is constructed from a variable 
of type tangent_vector. Upon construction, a geodesic is a function of type 
worldline, as defined in Section [5.4.2[ The code of the geodesic class is hsted 
in Section [A. 61 

The class geodesic maintains a list (cache) of all tangent_vectors found so 
far on the geodesic. The operator() routine of class geodesic, which takes tgnai 
as its only argument, uses the class bulirsch_stoer of Section H3] to attempt 
to numerically integrate fl6.3p from initial data in the cache. The particular 
initial data chosen is that whose affine parameter is nearest to tfinai- 



6.3 Parallel transport 

The operation of parallel transport represents the notion of transporting a 
vector along a curve while changing its direction as little as possible. It is 
defined in a similar way to a geodesicH 

A parallel transport is defined by a curve, and a tangent vector at a point 
on that curve. It then defines a unique tangent vector at each other point 
on the curve. On a chart, the components of the parallelly-transported 
tangent vector satisfy the equation 

where x'^{t) are the coordinates of the curve as a function of the curve pa- 
rameter t. 

Just as for geodesies, (16. 4p must in general be integrated on multiple 
charts to determine the tangent vector at a desired value t = tfinai of the 
curve parameter. The operation of parallel transport is represented in GR- 
workbench by the functor class paralleLtransport, which is constructed from 
a tangent_vector and a worldline. It is then a function with signature optional 
<tangent_vector> (double), representing the tangent vector as a function of 
the curve parameter t. The paralleLtransport class uses a similar algorithm 
to the geodesic class to integrate (16. 4p on any chart for which it is possi- 
ble, bisecting the interval of integration if integration cannot proceed on any 
chart. 



geodesic is, by definition, a curve whose tangent vector is the parallel transport of 
itself along the curve. 
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A parallelly-transported vector has a physical interpretation which makes 
it potentially useful in constructing numerical experiments: it is a fixed co- 
ordinate direction for a locally non-rotating physical observer who is moving 
on a geodesic. For locally non-rotating physical observers moving on non- 
geodesic world-lines, the operation with the corresponding physical interpre- 
tation is Fermi- Walker transport (see [T7], pages 47-49), which has not yet 
been implemented in GRworkbench. 



6.4 Implicitly-defined geodesies 

The methods of Section 16.21 allow the computation of the unique geodesic 
solving the initial value problem comprising fl6.3l) together with the initial 
coordinates a:* and the initial components of the tangent vector dx'^/dt. How- 
ever, there are ways other than the initial value problem to define a geodesic. 
Two physically important examples are discussed in this section. 

6.4.1 Unique connecting geodesies 

Around every point there is a neighbourhood such that, given two points 
within it, there will be a unique geodesic that intersects both points. The 
way to find this connecting geodesic is the topic of this section. The problem 
can be formulated in the following way: given two points a and b, find a 
tangent vector v & Ta such that the unique geodesic passing through a with 
tangent v also passes through b. If f is a solution to this problem, then 
so is av for any a 7^ 0; changing the value of a simply changes the affine 
parameter value at which the geodesic intersects b. 

The problem of finding the tangent vector v, up to scaling by a real 
number, can be thought of as determining which direction, in space and 
time, to launch a geodesic from a such that it 'hits' b. We solve this problem 
by minimising, over all possible directions at a, the amount by which the 
launched geodesic 'misses' b. To do this, we need a definition for the amount 
by which the geodesic misses — a real-valued function to minimise. 

The function / : — > M, which gives the amount by which the geodesic, 
launched from a with the given tangent vector, misses b, must satisfy certain 
properties. It must be zero for a geodesic which exactly intersects the point 
a, and strictly greater than zero otherwise; and it must be continuous, in the 
sense that, whenever a sequence of vectors f„ satisfy lim„^oo f{vn) = 0, then 
we must have f„ v, where v is an exact solution to the problem. 
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min_euclidean_separation 

A simple definition for the function /, satisfying the requirements listed 
above, is as follows: 

/(f ) = min min Sx^ = geodesic{v)(t)\c — b\c, (6.5) 

tgM charts C 

where geodesic(w) denotes the geodesic with tangent vector t> G Tq at a, \\ ■ \\ 
denotes the standard Euclidean norm on M", and we have used the notation 
that, for any point q and chart C, qc denotes the coordinates of q on C . That 
is, the distance between the curve geodesic (f) and the point h is defined as the 
closest they ever get, in the Euclidean norm, in the coordinates of any chart. 
The quantity is intended to be a small displacement in the coordinates 
of the chart C; in any case, it will certainly be zero if geodesic(f ) intersects 
h at affine parameter value t. 

The definition (16.51) is adequate, and was briefly employed in GRwork- 
bench, but it has a practical disadvantage: By using the Euclidean norm on 
R"', it effectively assigns equal importance to each of the coordinates. This is 
not ideal for some common coordinate systems. For example, consider, in the 
cylindrical coordinate system (t,r,(f),z) of (15. 5p . the point p = (0,10^^,0,0). 
Then the two points p+r = (0, 10"^+!, 0, 0) and = (0, 10'^, 1, 0) are equidis- 
tant from p in the sense of (16.51) . but p^r is much closer than p+0 to p in the 
sense of the standard fiat metric fl5.4p . essentially because the coefficient of 
the dr"^ term in fl5.4l) is 1, whereas the coefficient of the d(f)'^ term is = 10*. 

If the metric is diagonal, as above, then we can assign to each coordinate 

■2 

direction an approximate 'importance' equal to the coefficient of dx^ in 
the line element. If the metric is not diagonal then we diagonalise it at b, 
using the methods of Section I6.1.2[ and express the coordinate displacement 
6x^ in terms of the resulting orthonormal basis S of Tf,lf| 

6x'\b = B~^6x\ (6.6) 

Like the coordinate displacement 6x^\b € M" will depend on the chart C. 
The value is, in general, a better definition than for the amount 

by which geodesic(f )(t) 'misses' b\c, because it accounts for the difference in 
importance of the various coordinate directions at b. 
We rewrite f l6.5l) . using (16. 6p . as 

f{v) = min min \\6x'\b\\. (6.7) 

tgK charts C 



■^In (|6.6p . B is the matrix whose columns are the components of the orthonormal basis 
of Tf, on the chart C. 
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Definition fl6.7p is implemented in the routine min_euclidean_separation, which 
takes one argument of type worldline, and one argument of type point; it 
performs the minimisation of / over the curve parameter t using the one- 
dimensional minimisation routine brent_minimiser, of Section H75l 

Parameterising the search space 

We want to minimise the function f{v), fl6.7p . over the variable v E Ta- 
The tangent space Ta has dimension n, but, as already noted, f{av) = f{v) 
whenever a 7^ 0, and so the space to minimised over has dimension n — 1. 

In GRworkbench, the minimisation is performed in the following way: 
The vector v E Ta is expressed in terms of its components f * G M" in an 
orthonormal tangent basis B (Section 16.1.21) . Then, we minimise f{v) with 
f * ranging over the unit sphere in M", by parameterising the unit sphere 
by the n — 1 coordinates (6*1, ... , ^n-i) using the generalised spherical polar 
coordinate transformation, 

= sin 9i, 

m— 1 

= sin6'm Y\_ cos6'i, (1 < m < n), 

i=l 

n-l 

t;" = JJcos^i. (6.8) 

i=l 

The multi-dimensional minimisation of f{v) is thus performed over the n — 1 
variables (6*1, ... , 6'„_i). 

If the determined minimum value of f{v) is approximately equal to zero 
(in the sense of Section 14.1.11) . then the solution values (^1, . . . , 9n-i) of the 
minimisation problem define, via fl6.8p . the components of v in the or- 
thonormal tangent basis B, which in turn defines the solution vector v G T^, 
which finally defines, with a, the initial data for a geodesic intersecting a and 
b, as required. 

Implementation in GRworkbench 

The generalised spherical polar transformation fl6.8p is encoded in GRwork- 
bench in the routines from_polar and to_polar, both of which have signature 
nvector<double> (nvector<double>). The routine from_polar encodes (16. 8p . 
and to_polar encodes the inverse transformation to (16. 8p . The code for these 
routines is listed in Section IA.7I 
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The entire process of first solving the minimisation problem, 

min f(v), (6.9) 

by parameterising the space and minimising over the generalised spherical 
polar coordinates, and then constructing and returning the geodesic defined 
by the solution to (16.91) . is encapsulated in the routine connecting_geodesic 
of GRworkbench, which has signature optional<geodesic> (point, point). The 
optional mechanism is employed because it may not be possible to find the 
connecting geodesic; for example, the numerical minimisation of (16.91) may 
converge to a local, rather than a global, minimum, where f{v) ^ 0. The 
code of connecting_geodesic is hsted in Section 

The connecting_geodesic routine uses the routines to_polar and from_polar 
to perform the generalised spherical polar coordinate transformation, and 
the functor class powelLminimiser of Section 14.5.11 to perform the multi- 
dimensional minimisation. 

The minimisation class powelLminimiser requires an initial guess for the 
location of the minimum, around which it looks for an exact minimum; the 
guess supplied to powelLminimiser by connecting_geodesic is simply the coor- 
dinate difference between the two points a and h on some chart, transformed 
to the generalised spherical polar coordinates by the routine to_polar. This 
guess is good if the space-time curvature between a and h is small. 



6.4.2 Connecting null geodesies 

Given any world- line A(s) and a nearby point p, there will be two null 
geodesies which connect p with a point on A, corresponding to the intersec- 
tions of A with the the past and future null cones of p. These null geodesies 
are important because, if A is the world-line of a physical observer, they 
represent the world-lines of photons travelling to the event p from the ob- 
server, and from the event p to the observer. The determination of these null 
geodesies is the topic of this section. 

We solve the problem in a very similar way to the solution of the problem 
of Section l6.4.1l above: We minimise, over all null vectors v G Tp, the amount 
by which a geodesic launched from p with tangent vector v 'misses' the world- 
line A. There are two important differences between the two problems: We 
require a definition for the amount by which a curve misses another curve, 
analagous to the function / of (16. 7p which gives the amount by which a curve 
misses a point; and we only wish to minimise over null vectors in Tp, rather 
than all vectors in T„. 
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min_euclidean_separation of two curves 

We require a function g: Tp ^ M, analagous to / of (16.71) . which we can 
minimise to find the tangent vector at p of a null geodesic intersecting the 
point p and the world- line A. We define g in a similar way to /, as 



where the quantity Sx''\b is defined in terms of the quantity Sx^ as in (16. 6p . 
using an orthonormal tangent basis B at A(s), and 5x* is redefined as 



so that it is now a function of s, as well as t and C. 

We can summarise (I6.10p as follows: the distance between two curves is 
defined as the closest they ever get, in the Euclidean norm, in the coordi- 
nates of any chart. The definition (I6.10p is encoded in a specialisation of the 
GRworkbench routine min_euclidean_separation, which takes two arguments of 
type worldline, representing the space-time curves; it performs the minimisa- 
tion (I6.10p . over the two real parameters s and t, using the multi-dimensional 
minimisation routine powelLminimiser of Section [4.5.11 

Parameterisation of the null cone 

We want to minimise the function g{v), (I6.10p . over the variable v G Tp. The 
null subspace of Tp has dimension n — 1, and, as in Section [6.4.1l g{av) = g{y) 
whenever a 7^ 0, and so the space to be minimised over has dimension n — 2. 

The search space is parameterised in a similar way to that of Section [6.4.1t 
The vector f G Tp is expressed in terms of its components f * G M" in an 
orthonormal tangent basis B at p. Let the first vector in the basis B be the 
time-like eigenvector, and thus let the remaining eigenvectors be space-likejll 
The component thus represents the 'time-like part' of v, and the remaining 
n — 1 components v^, (3 = 2, . . . ,n, represent the 'space-like part' of v. Now, 
given any values for the v^, if we set 



then the vector v defined by the components f * is null, since the tangent basis 
B is orthonormal. Thus, to restrict our minimisation to the null space of Tp, 




(6.10) 




(6.11) 




(6.12) 



^In doing this, we implicitly assume that the metric is Lorcntzian, which is the usual 
case for physical applications of GRworkbench. 
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we minimise g over the components v^, and fix the remaining component f ^ 
using (Km . 

We minimise g{v) over the components G M"^^ by using the generahsed 
spherical polar coordinate transformation fl6.8p to obtain from the the 
coordinates (^i, . . . ,0n-2), which parameterise the unit sphere in ]R"~^, and 
then minimise g{v) over the n — 2 variables (^i, . . . , 6n-2)- 

As in Section [6.4.H if the minimum located value of g{v) is approximately 
equal to zero, then the solution values (6*1, ... , 9n-2) define the components 

via the generalised spherical polar coordinate transformation, and the v^, 
together with (16.121) . define the components of v in the tangent basis B, 
which in turn define the solution vector v G Tp, which finally, together with 
the point p, defines initial data for a solution geodesic intersecting both p 
and A. The geodesic is guaranteed to be null, due to fl6.12p . 



Implementation in GRworkbench 

The GRworkbench routine connecting_nulLgeodesic implements the process 
described above for minimising the function g{v) over all v in the null space of 
Tp, and constructing the resulting geodesic, using the to_polar and from_polar 
routines, and the powelLminimiser class. The signature of connecting_nulLgeodesic 
is optional<std::pair<double, geodesic>> (functional<optional<point> (double 
)>, point, double). The code of connecting_nulLgeodesic is hsted in Sec- 
tion |D 

The first and second arguments represent A and p, respectively. The 
third argument is an initial guess for the value of the parameter s of the 
world-line A, such that the connecting null geodesic will intersect A(s). This 
third argument is necessary for two reasons: There is otherwise no natural 
way for the connecting_nulLgeodesic routine to choose an initial guess for the 
values of the generalised spherical polar coordinates (^i, . . . , 6'„„2) to pass to 
powelLminimiser; and it permits a degree of control over which of the two 
possible connecting null geodesies (corresponding to either the backward or 
forward null cone of Tp) the connecting_nulLgeodesic routine will converge to. 

The return type of connecting_nulLgeodesic, optional<std::pair<double, 
geodesic>>, represents, in the first element of the std::pair, the parameter 
value s of the curve A at which the null geodesic intersects A; and in the 
second element of the std::pair, the null geodesic itself. By convention, the 
null geodesic returned by the routine connecting_nulLgeodesic intersects the 
curve A at the parameter value 1. 
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6.5 Conclusion 

Various tools useful for the simulation of physical situations have been im- 
plemented in GRworkbench. The tools are written within the functional 
framework of GRworkbench, allowing them to be easily interfaced with one- 
another to construct potentially complex physical models. Algorithms for 
the determination of implicitly-defined geodesies, in particular, demonstrate 
the numerical solution of an important physical problem using the numerical 
methods of Chapter [Hand the differential geometric framework of Chapter [51 



Chapter 7 

Analysis of a recent claim 



In this chapter we introduce and investigate a recent claim by Karim et al. [8] 
that the mass of the Milky Way can be determined using a small Michelson 
interferometer located on the surface of the Earth. After summarising their 
calculation in Section mi we analyse consequences of the physical model em- 
ployed by Karim et al. in Section [721 An alternative model, argued to be the 
correct one on physical grounds, is proposed and investigated in Section I7.3[ 
In Chapter [8] we describe numerical experiments performed in GRwork- 
bench using both models, and compare the results. 

7.1 Summary of the claim 

Employing a model metric of our galaxy, Karim et al. approximate the world- 
lines of the beam-splitter, end-mirrors, and connecting photons of an idealised 
Michelson interferometer located on the surface of the orbiting Earth. The 
proper time elapsed at the beam-splitter between the departure and return 
of photons along each interferometer arm is computed. 

The galaxy is modelled using a Kerr black hole metric. In Boyer-Lindquist 
coordinate^ {t, r, 9, (p), the Kerr metric takes the form 

ds"^ = gtt dt^ + 2gt^ dt dcp + dr"^ + gee d9^ + g^^ d(f)^. (7.1) 

The metric components gab depend on two parameters, m and a, which rep- 
resent, respectively, the mass and specific angular momentum]^ as measured 
from infinity, of the field source. Using the approximation employed by Karim 
et al, that a is small compared to m, and that m is small compared to the 

^See for example "7], page 161. 
^(angular momentum per unit mass) 
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components arr*^' 



radius of the orbit of the Earth about the centre of the galaxy, the metric 



a^sm^^-^ 
gtt = 5 -(1 - 2m/r), 

2mar sin^ 9 2m ^ n 
gtd> = 5 ^ a sm 9, 



gr 



p' 



gee = ^r"^ 



^ 1 — 2m/r ' 
+ 0^)2 - ^asin2 6' 



where 



P' 



/• 2 n i2 2 2 , -Z Z a 

4 = r — 2mr + a , p = r + a cos y. 



r^sin^^, (7.2) 



The world-hne of the beam-sphtter is modelled as a circular equatorial 
orbit about the centre of the galaxy: r = R, 9 = tt/2, and = 0o + ("^/-R)^; 
where R is the coordinate distance of the beam-splitter from the field centre, 
V is the coordinate speed of the beam-splitter, and v/R is the corresponding 
angular coordinate speed. The constant 0o is chosen to be zero. 

Karim et al. compute light travel times, to go up and back an interferom- 
eter arm, for three possible orientations of the interferometer arm: inward- 
radially directed, positive-0 directed, and positive-^ directed. Each arm is 
intended to have the same length L. 

The world-line of the end-mirror of the inward-radially directed arm 
(henceforth 'radial arm') is approximated as a circular orbit inside that of 
the beam-splitter: r = R — L, 9 = 7r/2, and = {v/R)t. The world-line of 
the end-mirror of the positive-0 directed arm (henceforth '0 arm') is approx- 
imated as a circular equatorial orbit which leads the beam-splitter in the 
direction by the angle $ = L/R: r = R, 9 = tc/2, and = $ {y/R)t. 
The world-line of the end-mirror of the positive-6' directed arm (henceforth 
'6' arm') is approximated as differing from that of the beam-splitter only in 
the 9 direction, again by the angle $: r = R, 9 = n/2 + <l>, = {y/R)t. 

The world-line of a photon travelling along an interferometer arm will in 
reality be a null geodesic which intersects the beam-splitter world-line, then 
intersects an end-mirror world-line, and finally intersects the beam-splitter 
world-line once again. To simplify the analytic calculation, Karim et al. 



■^Throughout, we use geometric units in which times are scaled by a factor c, and 
masses by a factor G/c^, so that physical quantities are measured in powers of metres. 
For example, angular momentum (kgm^ s~^) is measured in square metres. 
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make the approximation that the coordinates (r, 6, 0) are hnearly related 
along each photon world-line. The values of the remaining coordinate t for 
each world-line are fixed by requiring the tangent vector to the world-line to 
be null {ds = in (0))E1 

Explicitly, for photons travelling along the radial arm (where = 7r/2 
is constant by symmetry), d(j)/dr is assumed to be constant; for photons 
travelling along the arm, r and 6 are assumed to be constant; and for 
photons travelling along the 6 arm, r and d(p/dd are assumed to be constant. 

To summarise, Karim et al. make the following assumptions and approx- 
imations: 

1. Our galaxy is modelled by the Kerr black hole metric fl7.ip in the low 
angular- momentum approximation (17.21) . 

2. The Michelson interferometer is modelled in terms of the Boyer-Lindquist 
coordinates as described above. 

3. Photon world-lines are approximated as null curves in which the coor- 
dinates (r, 6', 0) are linearly related to one-another. 



7.1.1 Main results of the claim 

With the assumptions described above, Karim et al. solve for the coordinates 
of the arrival of a photon at the end-mirror, and for the return of the reflected 
photon to the beam-splitter. The t coordinate of the return event, scaled by 
the factor \/—gtt, gives the proper time elapsed at the beam-splitter. In 
terms of the dimensionless parameter /i = 2m /R and the coordinate speed 
V, Karim et al. find that the elapsed proper times for the radial, 0, and 6 
arms are, respectively. 



Tr = 2L 

T^ = 2L 
re = 2L 



. 1 5 2 1 2 
1 H — u u V + 

2^ 8^ 2 



1 9 1 9 



(7.3) 



The • ■ ■ denote terms of higher order in fi, v, and the parameter k = a/R. 



^The photon world-lines so defined, while null, will not, in general, be null geodesies. 
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field source 


2m (m) 


R (m) 


fi = 2m/R 


STre (s) 


Milky Way 

Sun 

Earth 


~ 1014 
~ 10^ 
~ 10-2 


~ 2.8 X 10^0 
~ 10" 
~ 6 X 10^ 


~ 10-6 
~ 10-s 
~ 10-s 


~ 6 X 10-15 
smaller 
smaller 



Table 7.1: Order of magnitude estimates of fl7.4p for various bodies with 
L = 10 cm, from [S]. 



For interferometry, the measurable quantity is the light travel time dif- 
ference between two arms. Karim et al. find that 



5TrB = Tr - Te ^ Ljl 

Sr^e = r^-Te^ 2Lv\ (7.4) 



and propose to determine /i (and hence the galactic mass m) by measuring 
the time differences (17. 4p . 

Karim et al. estimate the order of magnitude of the effect (17.41) due to 
the Earth, Sun, and Milky Way, for an interferometer of length 10 cm. The 
calculation is summarised in Table 17.11 The effect due to the Milky Way is 
found to be largest, with 

6Tre ~ 6 X 10-15 s. (7.5) 

Karim et al. conclude that the galactic mass can be determined by measuring 
5tj.q with a small interferometer. 



7.2 Theoretical analysis of the claim 

We now investigate properties of the physical model employed by Karim 
et al. In Section 17.31 we propose an alternative interferometer model, and 
investigate its properties. 

The main result in [8], upon which the proposed experiment depends, is 
the approximate light travel time difference (17.41) . It is independent of k and 
hence independent of a, the specific angular momentum of the gravitational 
field source. Thus, the result will be unchanged if the galaxy is instead mod- 
elled using a Schwarzschild black hole metric (setting a = in (17.21) ). In this 
case gtffi vanishes, and the algebra is simplified. We adopt this simpler model 
for the analytical calculations in Sections 17.2.11 and 17.31 and the numerical 
investigation of Chapter [HI 
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In discussing why the predicted time difference 6Tr0 is proportional to 
oc l/R, Karim et al. note that the proposed effect depends on the variation 



and suggest that it is therefore reasonable to expect an effect proportional 
to this potential. However, the variation of the potential over the volume of 
the interferometer will be approximately 

dm m 1 2 

Thus, it would seem that we should instead expect Sttq oc /i^. 



7.2.1 Properties of the coordinate-defined interferom- 



The interferometer of is defined in terms of the Boyer-Lindquist coordi- 
nates {t,r,6,(j)): The radial arm has coordinate length L in the r direction, 
and the 6 and arms have coordinate length ^ = L/R in the positive 6 and 
positive (j) directions, respectively. The justification for such a model is that, 
as R/2m oo, the metric components (17. 2p asymptote to those of the flat 
metric in spherical polar coordinates. 



and in that metric all of the arms of the coordinate-defined interferometer 
would have proper length L. 

Since geg and gcf,^ in (17.21) are equal to those in (17.71) . the 9 and arms 
of the coordinate-defined interferometer have proper length L. On the other 
hand, since Qrr in (17. 2p differs from that in (17.71) . the radial arm of the 
coordinate-defined interferometer does not have proper length L. In fact. 



^Karim et al. in fact describe 2rn/i? as the gravitational potential. In any case, since 
m/ R oc 2171/ R, the line of reasoning is unchanged. 




eter 




(7.7) 
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the proper length s of the radial arm is 



R-L 



r-R rR 
I R-L 

r-R / 

~ ' (l 

I R-L \ 

1 , R-L 
-2m In 



a/I — 2m /r 



dr 



12m\ , 

dr 

2 r J 



L 
L 



2 

1 L 
-2m — 

2 R 



R 



L + -L/i. 



(7.8) 



Consequences of model 

The proper length of the radial arm differs from L by an amount proportional 
to fi. The estimated time difference Stj.0 is also proportional to fi. This raises 
the possibility that the calculated value for 6Trg is due, at least in part, to the 
proper length difference between the r and 6 arms of the coordinate-defined 
interferometer. 

The total difference in proper length along and back each arm is 2(s— L) ~ 
Lfi. From fl7.4p . the lowest order term in 6Tre is also L/j,. This is exactly the 
time difference expected for an interferometer in fiat space, with arms of 
differing proper lengths s and L. We therefore conclude that the largest 
term in 6Trg, proportional to fi, is entirely due to the difference in proper 
lengths between the r and 6 arms of the coordinate-defined interferometer, 
and not to space-time curvature. 

Note that it does not follow from the above argument that there is no 
term proportional to ^ in the true physical value of 6Tre; it merely shows 
that, in the analysis of [8], the term proportional to /i is an artifact of the 
coordinate-dependent manner in which the interferometer is defined. Due to 
( 17.61) . however, we have good reason to believe that the lowest-order term in 
STrO is proportional to /x^, and not to fi. 



7.3 Geodesic-defined interferometer 

The problems resulting from the coordinate-dependent interferometer model 
of Karim et al. suggest that we should look for a cooidinaXe- independent 
model; we develop such a model in this section. Its properties are explored 
in Sections 17.3. 21 and l7. 3. 31 Along with the original model of Karim et al., this 
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alternative model is investigated numerically using GRworkbench in Chap- 
ter El 

7.3.1 Definition 

We begin by specifying the world-line of the beam-splitter in the same way 
as Karim et al: r = R, 6 = 7r/2, and = {v/R)t. Since this world-line 
will not, in general, be a geodesic|j it models an accelerating interferometer, 
rather than a freely-falling one. 

When deciding how to model the world-lines of the end-mirrors of each 
interferometer arm, the most obvious requirement is that the arms have 
length L. While the proper distance between two nearby points in a space- 
time may be defined as the proper length of the unique geodesic connecting 
them, it is a consequence of special relativity that there is no such observer- 
independent definition of the distance between two nearby world-lines. There 
is, however, a natural choice for a preferred observer: the beam-splitter, since 
proper time along the world-line of the beam-splitter is the physical quantity 
to be measured. 

With respect to a preferred observer, we can define the property of si- 
multaneity of two eventsJil Let b(r) be the world-line of the beam-splitter, 
where r is the proper time on b, let p be a point on b, let Tp be the tangent 
space of p, and let Aq G Tp be the tangent vector to b(r) at p. The vector 
Ao is the time direction of the beam-splitter at p. Let Sp be the space-like 
subspace of Tp orthogonal to Aq. The vectors in Sp are the space directions 
of the beam-splitter at p. An event q not on b is simultaneous with p if the 
unique geodesic connecting p and q has tangent v & Sp ai p. That is, q is si- 
multaneous to p if it is reachable from j9 by a (space-like) geodesic orthogonal 
to b. 

If q is simultaneous to p Eh then the distance between p and q is defined 
as the proper length of the space-like geodesic connecting them. 

Using the above definitions, we can construct an end-mirror world-line 
which is always a distance L from the beam-splitter. At each point p G b 
take a geodesic through p whose tangent vector Ai is orthogonal to Aq, and 
trace it out to proper length L. The end-point of this geodesic segment 
defines a point on the world-line of the end- mirror. To construct a second 
interferometer arm, take a second geodesic through p whose tangent vector 
A2 is orthogonal to both Aq and Ai, and trace it out to proper length L, 

^For each value of v there will be one value of R such that the world-line of the beam- 
splitter is, in fact, a circular equatorial geodesic. 

^For discussion regarding this definition of simultaneity see [1], pages 274-280. 
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defining the end-point as a point on the world-hne of the second end-mirror. 
A third interferometer arm can be similarly constructed by taking a third 
vector A3 which is orthogonal to Aq, Ai, and A2. 

Since the vectors (Ai, A2, A3) must be chosen for each p = b(r), they are 
functions of r. We require (Ai J\25 A3) to satisfy the following condition: each 
vector must vary continuousljo with r. This ensures that the orientation of 
the interferometer does not vary discontinuously. 

The vector Aq is fixed by our choice for the world-line b(r) of the beam- 
splitter, and 5'(Ao, Ao) = —1. We then choose Ai, A2, and A3 to model as 
closely as possible the same physical situation as Karim et al.u 

Xoocdt + ^d^, 
K 

Ai = —dr, 
A2 = de, 

X3 = d^ + giXo,d^)Xo. (7.9) 

To see that this set is orthogonal, observe that in the Kerr space-time 
{dt, dr, de) are mutually orthogonal, as are (c?<^, dr, de), while 

g{XQ, A3) = g{\o, d^ + g{\o, d^)\o) 

= g{Xo, d^) + g{Xo, d^)g{Xo, Aq) 
= 9iXo,d^) - g{Xo,d^) 

= 0. (7.10) 

Note that in the Schwarzschild space-time dt is orthogonal to d,/,, and so if 
V = then Aq is orthogonal to d^ and thus A3 = 9^. 

It remains to specify the world-lines of the photons connecting the beam- 
splitter to the end-mirrors. Now, from Section 16.4.21 given any world-line 
c and a nearby point p, there will be two null geodesies which connect p 
with a point on c, corresponding to the intersections of c with the the past 
and future null cones of p. Thus, for each interferometer arm, we let the 
world-line of the outgoing photon be the (locally unique) future directed null 
geodesic joining the origin event O to some point q on the world-line of the 
end-mirror, and we let the world-line of the returning photon be the future 
directed null geodesic joining q to some point r on b. 

For each arm, the proper length of b between the origin event O and 
the return event r is the time experienced by the beam-splitter between 
the departure and return of a photon travelling along that arm. The point 

®The components of each vector must be continuous in any coordinate system. 
^The shorthand notation d^^i represents the coordinate basis vector d/dx^ . 
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r will in general be different for each interferometer arm, and the proper 
length along b between two such points gives the measurable light travel 
time difference between the corresponding interferometer arms: 6Tre, 6Tr,j,, or 

7.3.2 Comparison with the coordinate-defined inter- 
ferometer 

The geodesic-defined interferometer has the following properties: 

1. The arms are of proper length L. 

2. The arms are straight, in the sense of a geodesic being the straightest 
possible line in a curved space. 

3. At their point of intersection, the arms are orthogonal to: 

(a) one-another; 

(b) the world-line of the beam-splitter. 

We have seen in Sections 17.11 and 17.21 that properties [1] and [2] are not shared 
by the coordinate-defined interferometer of Karim et al. . 

Property [3a| is shared by the coordinate-defined interferometer, because 
the tangent vectors to the arms are 9^, de, and S^, which are mutually or- 
thogonal. Property [3b] does not hold in general because, when f 7^ 0, the 
tangent vector to the world-line of the beam-splitter (equal to Aq, above) is 
not orthogonal to d^] and because in the Kerr space-time dt is not orthog- 
onal to dif). In the special case of the Schwarzschild space-time with w = 0, 
property [3b] does hold for the coordinate-defined interferometer. 

7.3.3 Estimate of fight travel time 

In this section we estimate Tr for the geodesic-defined interferometer, for the 
simplest case of f = in the Schwarzschild space-time. We will find that the 
result differs from 2L by an amount proportional to /x^, in contrast to the 
corresponding result (17.31) for the coordinate-defined interferometer. 

From symmetry it follows that the world-line of the outgoing radial light 
ray has constant 6* = 7r/2 and constant (j). Since the world-line is null, (is = 
along it. Thus, from (17.11) and (17.21) . with a = 0, 



= gtt dt^ + Qrr dr"^, 



(7.11) 
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where 

Let the time coordinate of the photon leaving the beam-sphtter be t = 0, 
and let the time coordinate of the photon reflecting at the mirror he t = t,-. 
Then, since the space-time is static and time-reversible the time coordinate 
of the return of the photon to the beam-splitter is 



2U, (7.12) 



in terms of which 
From (17. lip we have 



Tr = 2t,v/^. (7.13) 



R 



^dr, (7.14) 
i?-A V 9tt 

where A is the coordinate distance on the r axis corresponding to a proper 
length L. 

Relation between coordinate length and proper length 

To find an expression for A in terms of L, we first find L in terms of A: 

R rR I 



y/g^dr= / — dr. (7-15) 

R-A Jr-a V 1 - 2m/r 

The solution to this integral can be expressed in closed form, but we only 
require the first few terms in A^, = A/2m. Using Mathematica we obtaiio 



R,-l 4v^(i?, - 1)3/2 * 

+ ^nr^ + 7-16 

24i?P(i?, - 1)5/2 



where R^ = R/2m = l//i and L^, = L/2m. We can inverto this series to 
obtain a series for A^, in terms of L*. We begin by rewriting fl7.16p as 

L^ = aiA, + a2Al + a3Al + --- , (7.17) 



^"^The Schwarzschild space-time is static because the metric is independent of t, and 
time-reversible because it is invariant under the exchange t — s- —t, dt — > —dt. The Kerr 
space-time is thus static but not time-reversible. 

^^Mathematica input: Simplify[Series[lntegrate[l / Sqrt[l — 1 / r], {r, R — Delta, R}], 
{Delta, 0, 3}]] 

^^The general process of finding a series which is the inverse function of another series 
is called series inversion or series reversion. 
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and then writing a general series for A* in terms of L^: 

A,=hL, + b2Ll + hLl + --- . (7.18) 

Substituting (17.171) into (I7.18p . equating powers of A^,, and solving for the bi 
yieldfl 

&2 -- 



a 



3 ' 

1 



63 = MZ^. (7.19) 



The series for A* in terms of is thus 



-R* — 1^ -'-7-2 1 r^* ^ ^ t3 



Solution 

We evaluate the integral (I7.14p for tj.* = tj-/2m using the reduced variable 
r^, = r/2m: 

tr*= [ ^T^dn = A, + ln „ \ ■ (7.21) 



,-A, 1 — 1/?"* i?* — 1 — A* 

Substituting (17.201) and (17.211) into (17.131) and expanding in powers of and 
yields, after some simplification, 

Ll LI/3-LI/4 
2m~ * 2Ry Rf 

= 2L,-^L^ + ... , (7.22) 

or 

= 2L - + ■ ■ ■ . (7.23) 

Thus, Tj. differs from 2L by an amount proportional to /i^, in agreement 
with the argument of (17. 6p . (17.230 . along with (17. 3p . will also be useful in 
validating the numerical analysis of Chapter [HI 



^^For a formula for the general coefficient 6„, see, for example, [TTj, page 412. 
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The quantity mL"^ / in fl7.23p differs from the corresponding quantity 
2mL/R from the analysis of Karim et al. (17.31) by a factor of L/2R. For 
L = 1 m and R = 8kpc0 we have L/2R ^ 2 x 10^^^. Thus we might expect 
a change in the time difference estimate flT.Sp of roughly a factor of 10~^^, 
so that 5Tre ~ 10~^^ s, which is too small to detect with current methods. 
An accurate estimate of the time difference ST^e for the geodesic-defined 
interferometer is obtained numerically in Chapter [HI 



7.4 Intermission 

Because it is defined explicitly and simply in terms of the Boyer-Lindquist 
coordinates {t,r,6,(j)), the coordinate-defined interferometer of Karim et al. 
is more susceptible to analytic methods than the geodesic-defined interfer- 
ometer of Section 17.31 Nonetheless, to keep the algebra manageable, various 
approximations were necessarily employed in [S]. In particular, by approx- 
imating null geodesies as null curves in which the coordinates (r, 6, 0) are 
linearly related, Karim et al. completely avoid the geodesic equation in their 
analysis!^ 

The geodesic-defined interferometer, on the other hand, is defined explic- 
itly terms of space-like geodesies, and so an analysis of it akin to that of 
[H] would be even more complicated. We do, however, have the methods of 
Chapter [6] at our disposal. In Chapter [8] we directly simulate both interfer- 
ometers, bypassing the algebraic complexities of the metric and the geodesic 
equation. By performing a range of numerical experiments, we can charac- 
terise the behaviour of both interferometers in terms of the parameters R, 
L, and v. 



^''The estimate for R is taken from [3], page 917. 
Similarly, the analysis of Section 17.3.31 was relatively simple because the radial 
geodesies were easily found via the symmetries present in the special case v = a = 0. 



Chapter 8 

Numerical investigation of the 
claim 

Using the methods of Chapter [6], the coordinate-defined interferometer of 
Karim et al. and the geodesic-defined interferometer of Section [7^ were sim- 
ulated in GRworkbench, in the Schwarzschild space-time. In Section 18.11 the 
modelhng of the interferometers in GRworkbench is described. The results 
of the numerical experiments are presented in Section 18.31 In Section 18.41 
the results for the geodesic-defined interferometer are used to obtain a new 
estimate for the size of the predicted effect on Earth due to the Milky Way. 
Conclusions are drawn in Section 18. 5[ 

The motivation for the experiments was twofold: Under the assumption 
that the geodesic-defined interferometer is more physically realistic than the 
coordinate-defined interferometer of [8], we aimed to obtain a new estimate 
on the size of the effect Sr^e, in order to determine whether the Milky Way 
can in fact be weighed with a small interferometer on Earth; and we aimed to 
verify the analysis of the coordinate-defined interferometer made in [8]. By 
directly simulating the coordinate-defined interferometer, we can bypass the 
approximations necessary in an analytic argument, including the approxima- 
tion of light rays as certain (non-geodesic) null curves, and thus determine 
the extent to which those approximations affected the final result of Karim 
et al. 

8.1 Modelling the interferometers 

In this section we describe how to simulate the two interferometer mod- 
els defined in Chapter [71 using the tools for numerical experimentation de- 
scribed in Chapter El The coordinate-defined interferometer of Karim et 
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al. is constructed in terms of straight lines in coordinate space, using the 
coordinate_line tool of Section [6.1.3l while the geodesic-defined interferometer 
also makes use of the geodesic functor class of Section 16.21 For both inter- 
ferometers, null geodesies, representing photon world-lines, are determined 
using the implicit methods of Section 16.41 

Each interferometer model depends on the three parameters R, L, and v, 
corresponding, respectively, to the coordinate distance of the beam-splitter 
from the field centre, the interferometer arm length, and the coordinate speed 
of the beam-splitter (Section 17.11) . The important physical quantities ob- 
tained from each simulation are the light travel time differences STro, STr^, 
and St0^, which are arc lengths along the world-line of the beam-splitter. By 
simulating each interferometer model for a wide range of values of R, L, and 
V, the effect of each parameter on the light travel time differences can be 
characterised. 

8.1.1 Beam-splitter world-line 

For both interferometer models, the world-line of the beam-splitter is mod- 
elled as a circular equatorial orbit, which is a straight line in the Boyer- 
Lindquist coordinates (t,r,9,(f)). The world-line satisfies (Section |7. II) 

t = s, r = R, e = 7i/2, (j) = (j)o+ {v/R)s, (8.1) 

where s is a curve parameter; the tangent vector to the curve (18. ip every- 
where has the components (1, 0,0,v/R). However, the parameter s does not 
correspond to the proper time r of the beam-splitter, because the vector Aq 
with components Aq = (l,0,0,f/i?) does not satisfy metric(Ao, Aq) = — 1. 
We normalise Aq using the routine normalise of Section [6. l.lj and use the re- 
sulting vector u to construct a coordinate_line whose parameter is the proper 
time r. The arbitrary constant 0o is chosen to be 7r/2. Note that Aq as 
defined here is simply the Aq of (17. 9p . 

The following code fragment demonstrates the construction of the beam- 
splitter world-line in GRworkbench: 

// construct the point representing the origin event 

nvector<double> origin_coordinates = make_vector(0, R, half_pi, half_pi); 
point origin(a, c, origin_coordinates); 

// construct the world— line of the beam— splitter 

nvector<double> coordinate_direction = make_vector(l, 0, 0, v / R); 

tangent_vector beam_splitter_tangent = normalise(tangent_vector(origin, c, 

coordinate_direction)); 
worldline beam_splitter_worldline = coordinate_line(beam_splitter_tangent, c); 
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The variable c is assumed to be of type chart, representing a chart which 
uses the Boyer-Lindquist coordinates, and the variable a is assumed to be 
of type atlas, representing the Schwarzschild space-time encoded in GRwork- 
bench. After execution of the code fragment, above, the beam-splitter world- 
line, represented by a function of type worldline (Section I5.4.2I) . is stored 
in the variable beam_splitter_worldline, and the argument to the function 
beam_splitter_worldline, of type double, corresponds to the proper time of 
the beam-splitter. 

Note that the coordinate_line on the last line of the code fragment, above, 
is constructed from a tangent_vector and a chart; the information regarding 
the origin point is contained in the context routine of the tangent_vector class; 
see Section I5.4.3I 

8.1.2 End-mirror world-lines 

Both interferometer models have all parts of the interferometer orbiting the 
field centre at a constant value of the r coordinate. Hence, the end-mirror 
world-lines, like the beam-splitter world-line, have tangent vectors whose 
components are proportional to (1, 0, 0, f /-R). The only difference between 
the construction of an end-mirror world-line in GRworkbench, and the con- 
struction of the beam-splitter world-line in the code listing, above, will be 
the definition of the variable origin of type point. 

Coordinate-defined interferometer 

For the coordinate-defined interferometer, the origin events of the end-mirrors 
are defined simply in terms of the Boyer-Lindquist coordinates. For the 
inward-radial arm, the origin event has coordinates (0, i? — L, 7r/2, 7r/2); for 
the positive-(/) arm, the origin event has coordinates (0, i?, 7r/2, 7r/2 + L/R); 
and for the positive-^ arm, the origin event has coordinates {0,R,tt/2 + 
L/R,7r/2). The following code fragment demonstrates the construction of 
the end-mirror world-lines in GRworkbench: 

// (choose one of the following three lines) 

nvector<double> mirror_origin_coordinates = make_vector(0, R — L, half_pi, 

half_pi); // inward— radial arm 

nvector<double> mirror_origin_coordinates = make_vector(0, R, half_pi + L / 

R, half_pi); // positive— theta arm 

nvector<double> mirror_origin_coordinates = make_vector(0, R, half_pi, 

half_pi + L / R); // positive— phi arm 



// construct the point representing the origin event 
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point mirror_origin(a, c, mirror_origin_coordinates); 
// construct the world— line of the end— mirror 

nvector<double> coordinate_direction = make_vector(l, 0, 0, v / R); 
tangent_vector mirror_tangent = normalise(tangent_vector(mirror_origin, c, 

coordinate_direction)); 
worldline mirror_worldline = coordinateJine(mirror_tangent, c); 

The only significant difference between this code fragment, and the code frag- 
ment demonstrating the construction of the beam-splitter world-line, above, 
is in the definition of the coordinates of the origin point. 



Geodesic-defined interferometer 

As described in Section I7.3.H the origin events for the end-mirrors of the 
geodesic-defined interferometer are the end-points of space-like geodesies of 
length L emanating from the origin event of the beam-splitter, and orthogo- 
nal to the world-line of the beam-splitter. The tangent vectors of the space- 
like geodesies at the origin event are the mutually orthogonal vectors Ai, A2, 
and A3, of ra . 

The vectors Ai, A2, and A3 are obtained from the coordinate basis vectors 
dr, dg, and d^f, by using the orthonormalise routine of Section [B.l.li Specifi- 
cally, Ai is defined as the orthonormalisation of dr with respect to the tangent 
Ao to the world-line of the beam-splitter; A2 is defined as the orthonormal- 
isation of dg with respect to both Aq and Ai (obtained by two applications 
of orthonormalise); and A3 is defined as the orthonormalisation of c?^ with 
respect to Aq, Ai, and A2. This process is equivalent to applying the Gram- 
Schmidt process (see for example [S], page 399) to the vectors Aq, dr, dg, and 
d^, to obtain an orthonormal basis for the tangent space at the origin. 

The following code fragment demonstrates the construction of the end- 
mirror world-lines of the geodesic-defined interferometer in GRworkbench: 

// coordinate basis vectors 

tangent_vector r (mirror_origin, c, make_vector(0., —1., 0., 0.)); 
tangent_vector theta (mirror_origin, c, make_vector(0., 0., 1., 0.)); 
tangent_vector phi (mirror_origin, c, make_vector(0., 0., 0., 1.)); 

// gram— Schmidt process 

tangent_vector radiaLmirror_direction = orthonormalise( 

r, beam_splitter_tangent); 
tangent_vector theta_mirror_direction = orthonormalise(orthonormalise( 

theta, beam_splitter_tangent), radiaLmirror_direction); 
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tangent_vector phi_mirror_direction = orthonormalise(orthonormalise( 
orthonormalise( 
phi, beam_splitter_tangent), 
radiaLmirror.direction), theta_mirror_direction); 

// construct the space— like geodesic representing the interferometer arm 

// (choose one of the following three lines) 

worldline interferometer_arm = geodesic(r_mirror_direction); 

worldline interferometer_arm = geodesic(theta_mirror_direction); 

worldline interferometer_arm = geodesic(phi_mirror_direction); 

// determine the point representing the origin event of the end— mirror 
point mirror_origin = interferometer_arm(L); 

/ / construct the world— line of the end— mirror 

nvector<double> coordinate_direction = make_vector(l, 0, 0, v / R); 
tangent_vector mirror_tangent = normalise(tangent_vector(mirror_origin, c, 

coordinate_direction)); 
worldline mirror_worldline = coordinateJine(mirror_tangent, c); 

The difference between this code fragment, and the corresponding code frag- 
ment for the construction of the coordinate-defined interferometer, is in the 
definition of the origin event for the end- mirror — the variable mirror_origin. 
For the geodesic-defined interferometer, above, it is constructed in terms of a 
space-hke geodesic from the mirror_origin event, whereas, for the coordinate- 
defined interfermeter, it was constructed exphcitly in terms of the Boyer- 
Lindquist coordinates. 

8.1.3 Photon world-lines 

In Sections 18.1.11 and 18.1.21 the origin event mirror_origin, from which pho- 
tons are emitted, and the end-mirror world-hnes (mirror_worldline in the code 
fragment above), with which the photons must intersect, were defined. This 
is sufficient information for the apphcation of the method of Section 16.4.21 to 
obtain null geodesies representing the world-lines of outgoing photons. 

Once the outgoing geodesies have been obtained, their points of intersec- 
tion with the end-mirror world-lines define reflection events. The reflection 
events, together with the beam-splitter world-line, beam_splitter_worldline, 
constitute sufficient information to again apply the method of Section 16.4.21 
to obtain null geodesies representing the world-lines of ingoing photons. 

The points of intersection of the ingoing geodesies with the world-line of 
the beam-splitter will occur at various values of the world-line parameter r. 
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the proper time of the beam-sphtter. The difference between these values of 
r define the hght travel time differences 6Tre, Srr^, and 6t0^, which are the 
quantities to be obtained. 

The following code fragment demonstrates the application of the routine 
connecting_nulLgeodesic of Section [6.4.21 to determine the light travel time for 
one interferometer arm: 

geodesic outward_ray = connecting_nulLgeodesic(beam_splitter_origin, 

mirror_worldline, L)— >second; 
point reflection = outward_ray(l); 

double light_traveLtime = connecting_nulLgeodesic(reflection, 
beam_splitter_worldline, 2 * L)— >first; 

In the first line, the routine second obtains the second element of the std:: 
pair<double, geodesio returned by the routine connecting_nulLgeodesic (see 
the end of Section I6.4.2p . In the second line, we make use of the con- 
vention that the null geodesic returned by connecting_nulLgeodesic inter- 
sects mirror_worldline at parameter value 1. In the third line, the routine 
first obtains the first element of the std::pair<double, geodesio returned by 
connecting_nulLgeodesic, which corresponds to the parameter r of the world- 
line of the beam-splitter at which the ingoing photon arrives. 

Note that the third argument to connecting_nulLgeodesic, an initial guess 
for the parameter value of the curve at which the null geodesic will intersect, 
is chosen to be L for the outgoing ray intersecting with the end-mirror world- 
line, and 2L for the ingoing ray intersecting with the beam-splitter world- 
line. These guesses correspond to the exact points of intersection for an 
interferometer in flat space, where the light travel time will be L to reach 
the mirror, and 2L to return to the beam-splitter; they are good guesses if 
the space-time curvature is small in the region of interest. 

Figure 18.11 shows the coordinate-defined interferometer modelled in GR- 
workbench, as described in this section. There are 5 interferometer arms: 
inward-radial, outward-radial, positive-0, negative-0, and positive-^. (By 
symmetry, the negative 6 arm has the same light travel time as the positive 
6 arm.) The photon world-lines, determined by connecting_nulLgeodesic, are 
visible for both of the radial arms and the positive-6' arm. 

8.2 Experiment 

Using the methods of Section 18. ![ we can simulate either the coordinate- 
defined interferometer of Karim et ai, or the geodesic-defined interferometer 
of Section 17731 for any values of the parameters R, L, and v. Because physical 
values of L/R are smaller than 10~^^, the precision of the double type in 
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Figure 8.1: The coordinate-defined interferometer witli 5 ortliogonal arms, 
simulated in GRworkbench. 
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C++, it is not possible to directly simulate an interferometer on Earth under 
the influence of the galactic gravitational field. However, by simulating the 
interferometer for a wide range of values of R, L, and v, the dependence of 
the light travel time difference on each parameter can be discovered, and the 
effect at Earth due to the galactic gravitational field can be predicted. 

Appendix [B] lists the code of the numerical experiment performed in 
GRvjorkbench to characterise each of the interferometer models. The sim- 
ulation of the coordinate-defined interferometer is represented by the class 
karimJnterferometer, and the simulation of the geodesic-defined interferom- 
eter is represented by the class geodesic_interferometer. The reflect routine 
of each class performs the simulation of the corresponding interferometer; 
it takes three arguments of type double, representing the values of the di- 
mensionless parameters R^ = R/2m, = L/2m, and v, where 2m is the 
Schwarzschild radius for a black hole of mass m. 

The reflect routine computes the light travel times r^, Tg, and r^, as 
described in Section 18.11 and takes their difference to form the travel time 
differences Srrg, St„i,, and ^r^^. The computed travel time differences are in 
units of 2m. 

For each interferometer, 5 experiments were performed, with each ex- 
periment comprising many calls to reflect, that is, many simulations of the 
interferometer. The 5 experiments were 

1. v = 0, L^ = 1,3<R^< 50, and 

2. v = 10-2, = 1, 3 < i?, < 50 (varying i?,); 

3. v = 0, R^ = 10, 10^2 < < 6, and 

A. v = IQ-^, R^ = 10, IQ-^ <L^<6 (varying L,); and 
5. R^ = 10, = 1, 10-3 <v <0.5 (varying v). 

In the experiments, i?* was varied over 17 values in a geometric progression 
starting with R^ = 3; L^, was varied over 37 values in a geometric progression 
starting with L^, = 10"^; and v was varied over 37 values in a geometric 
progression starting with v = 10~^. Thus, each interferometer model was 
simulated for a total of 145 different sets of values for the parameters i?*, L^, 
and V. 

8.3 Results 

In this section we present the results of the numerical experiments described 
in Section I8l2l 
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8.3.1 Validation 

An analytic calculation for the light travel time along the radial arm of the 
geodesic-defined interferometer, for the special case v = 0, was made in 
Section 17.3.31 resulting in a power series expansion in L^, and for the 
travel time r^, ( \7.22\\ . This travel time was compared with the values for 
Tr obtained in the numerical experiments of Section 18.21 for various values 
of i?* and L^. In all cases the numerical experiment results were found to 
agree with the analytic calculation in the first 8 or 9 significant figures. The 
relative precision used by the approx_equal mechanism of Section [4.1.11 was 
10~^ for the numerical experiments described in this chapter. 

The case f = is not special from the point of view of the numerical 
differential geometric engine of GRworkbench. It can thus be extrapolated 
that the light travel times determined by the numerical methods of this 
chapter when v ^ are also as accurate as permitted by the relative precision 
of the numerical methods. 

8.3.2 Varying orbital radius 

Figures 18.21 and 18.31 show the light travel time differences for Experiment [T], 
of Section 18.21 for the coordinate-defined interferometer and the geodesic- 
defined interferometer, respectively. Note the logarithmic axes on these plots, 
and all plots in this section. 

In all figures in this section, three sets of data are plotted, corresponding 
to the light travel time differences between the three pairs of interferometer 
arms: r-6, r-(j), and 6-(j). 

The data for the r-6 time difference coincides with the data for the r- 
(f) time difference on Figures 18.21 and 18.31 because, when v = 0, the (p aiid 
6 arms are equivalent, owing the spherical symmetry of the Schwarzschild 
space-time. 

The relative precision of the numerically determined light travel time 
differences is at best 10~^; we see from Figures and that the 6-(f) time 
differences are well below the numerical precision limit — they are effectively 
zero. This is to be expected because, since the 9 and cf) arms are equivalent 
when V = 0, the light travel time along them should be exactly the same 
(within the numerical precision). 

From Figure [H^ for large values of R^^, the slope of the r-6 time difference 
data is very close to —1 on the logarithmic scale, corresponding to the travel 
time difference 6Tre being proportional to 1/i?* for the coordinate-defined in- 
terferometer. This 1/i?* scaling is in agreement with the calculation (17.41) of 
Karim et al. and, comparing the values of the r-9 data in Figure [8l2] with the 
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Figure 8.2: Light travel time difference for the coordinate-defined interfer- 
ometer, for various values of R^, with fixed L^ — 1 and v — 0. 
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Figure 8.3: Light travel time difference for the geodesic-defined interferome- 
ter, for various values of R^, with fixed = 1 and v = 0. 
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Figure 8.4: Light travel time difference for tlie coordinate-defined interfer- 
ometer, for various values of i?*, witli fixed L^, = 1 and v = 10^^. 



predicted travel time differences, fl7.4p is found to be accurate to several sig- 
nificant figures. Thus, the analysis of the coordinate-defined interferometer 
by Karim et al. is validated. 

For large values of R^,, the slope of the r-9 time difference data for the 
geodesic-defined interferometer (Figure 18.31) is found to be very close to —2 
on the logarithmic scale, corresponding to the travel time difference ^r^e 
being proportional to 1/Rl- This is in agreement with the argument (17. 6p of 
Section 17. 2[ 

Figures 18.41 and 18.51 show the light travel time differences for Experi- 
ment [21 for the coordinate-defined interferometer and the geodesic-defined 
interferometer, respectively. The physical situation modelled in producing 
these plots differs from that of Figures [F!2] and [8731 only in the interferometer 
coordinate speed v being non-zero for these plots. 

The r-6 data and the r-(f) data of Figures 18.41 and 18.51 do not differ sig- 
nificantly from the corresponding data for Figures 18.21 and 18.31 despite the 
non-zero interferometer coordinate speed. In particular, the data for the r-6 
time differences still coincides with data for the r-0 time differences, despite 
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Figure 8.5: Light travel time difference for the geodesic-defined interferome- 
ter, for various values of i?*, with fixed = 1 and v = 10~^. 
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the two arms 6 and being no longer equivalent. The coincidence of these 
two data sets is, in fact, a feature of all the plots in this section. 

The 6-(j) time difference data for the coordinate-defined interferometer 
(Figure [8. 41) can be seen to be roughly independent of i?*, for large values 
of -R^,. This is in agreement with the estimate (17.41) of the time difference 
6Tgfp of Karim et ai. Once again, examining the data comprising Figure [83, 
it is found to be in agreement with the estimate (17.41) in the first several 
significant figures, validating the analysis of Karim et al. 

Interestingly, for large values of R^,, the 6-(f) time difference data for the 
geodesic-defined interferometer has a slope very close to —3 on the loga- 
rithmic scale, corresponding to the time difference 5re(0 being proportional 
to 1/Rl- Thus, while the 6-(f) time difference is already smaller than the 
r-6 time difference on Figure 18.51 by several orders of magnitude, at physi- 
cal values of i?* (i?* > 10^), it will be comparatively even smaller. This is 
in contrast with the situation for the coordinate-defined interferometer: On 
Figure [8^ it would appear that, if we extrapolate the data to physical values 
of i?*, we might enter a regime where the 6-(j) time difference is larger than 
the r-6 time difference. 

8.3.3 Varying interferometer length 

Figures 18.61 and 18.71 show the light travel time differences for Experiment 0, 
for the coordinate-defined interferometer and the geodesic-defined interfer- 
ometer, respectively. 

As with the other experiment with v = (Experiment [T]), and as ex- 
pected, the 6-(f) time difference data is everywhere zero, within the numerical 
precision. 

For small values of L^,, the r-6 data for the coordinate-defined inter- 
ferometer has slope very close to 1 on the logarithmic scale of Figure 18.61 
corresponding to the travel time difference ^r^.^ being proportional to L*. 
Again, the scaling is in agreement with the estimate (17.41) of Karim et al. 

For the geodesic-defined interferometer, for small values of L*, the r-6 
data has slope very close to 2 on the logarithmic scale of Figure 18.61 corre- 
sponding to the travel time difference dT^e being proportional to L^. 

Figures 18.81 and 18.91 show the light travel time differences for Experi- 
ment m for the coordinate-defined interferometer and the geodesic-defined 
interferometer, respectively. 

As with Experiments [Hand [21 there is no significant difference between the 
r-d data of Figures 18.81 and 18.91 and the corresponding data from Figures 18.61 
and[H21 
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Figure 8.6: Light travel time difference for the coordinate-defined interfer- 
ometer, for various values of L*, with fixed = 10 and v = 0. 
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Figure 8.7: Light travel time difference for the geodesic-defined interferome- 
ter, for various values of L*, with fixed R^ = 10 and v = 0. 
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Figure 8.8: Light travel time difference for the coordinate-defined interfer- 
ometer, for various values of L*, with fixed — 10 and v — 10~^. 
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Figure 8.9: Light travel time difference for the geodesic-defined interferome- 
ter, for various values of L*, with fixed R^ = 10 and v = 10~^. 
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Figure 8.10: Light travel time difference for the coordinate-defined interfer- 
ometer, for various values of v, with fixed i?* = 10 and = 1. 

For small values of L^,, the slope of the 6-(j) data on the logarithmic scale 
of Figure 18.81 is very close to 1, corresponding to the travel time difference 
5r6)(^ being proportional to for the coordinate-defined interferometer. This 
scaling is in agreement with the calculation (17.41) of Karim et al. 

Almost all of the 6-(f) data for the geodesic-defined interferometer (Fig- 
ure 18. 9p are near or below the relative precision of the numerical methods, 
10~^, and so no reliable conclusions can be drawn about it. Based on the few 
reliable data points, which are unfortunately at large (non-physical) values 
of L^, we might conjecture an Ll dependence of 6Tg^ on L^,, consistent with 
the scaling of Sr^e, since the slope of the valid 9-(j) data points is roughly 2. 



8.3.4 Varying interferometer coordinate speed 

Figures [8.101 and 18.111 show the light travel time differences for Experiment [5l 
for the coordinate-defined interferometer and the geodesic-defined interfer- 
ometer, respectively. 

The most important property of these plots is that, for both interferome- 
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Figure 8.11: Light travel time difference for the geodesic-defined interferom- 
eter, for various values of v, with fixed i?* = 10 and = 1. 
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ter models, for small values of v, the r-6 time difference data are independent 
of V. For the coordinate-defined interferometer, this result is in agreement 
with the estimate fl7.4p of Karim et ai. For the geodesic-defined interferom- 
eter we conclude that, for physical values of w (f ~ 10~^), the travel time 
difference 6Trg is independent of v. 

For small values of v, the slope of the 6-(f) data on the logarithmic scale 
of Figure [5.101 is very close to 2, corresponding to the travel time difference 
6Tg^ being proportional to f ^ for the coordinate-defined interferometer. This 
scaling is in agreement with the calculation (17.41) of Karim et al. . 

The slope of the O-cj) data for the geodesic-defined interferometer (Fig- 
ure 18.111) is also very close to 2 for small values of although it should be 
noted that the first few data points are near or below the relative precision 
10~^ of the numerical methods employed. 

The unusual behaviour of the r-9 data on Figures [8. 101 and 18. lll for values 
of V approaching unity is simply due to the light travel time difference passing 
through zero on the logarithmic axes. Because f is a coordinate speed, if it 
is increased beyond approximately unity, then the world-lines of the various 
parts of the interferometer will become space-like, which is certainly not 
physical. 

8.3.5 Summary 

The results of the all the numerical experiments simulating the coordinate- 
defined interferometer were in agreement with the estimated light travel time 
differences (17. 4p of Karim et al. Thus, the analysis of the coordinate-defined 
interferometer in [S] was validated. 

The light travel time differences for the geodesic-defined interferometer 
were investigated as a function of the dimensionless parameters i?*, L*, and 
V. The largest travel time difference was Sr^e (or ^r^,/,), which was found 
to be proportional to Ll/Rl, independent of v, for small values of L^, large 
values of i?*, and small values of v. 

8.4 Estimate of physical effect 

In this section we employ the relation 

6Tre cc Ll/Rl (8.2) 

for the geodesic-defined interferometer, which was discovered by numerical 
experimentation in Section 18.31 to estimate the size of the light travel time 
difference 5T,.e for a 1 m interferometer on Earth. Analagous to Table 17.11 we 
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field source 


M(kg) 


R (m) 


5tsi (s) 


Earth 
Sun 

Milky Way 


5.97 X 10^4 
1.99 X lO^o 
2 X 10^^ 


6.38 X 10^ 
1.50 X lOH 
2.5 X 10^0 


3.5 X 10-25 
2.09 X 10-28 
8 X 10-^^ 



Table 8.1: Estimates of dr^i for various bodies with L = 1 m, for the geodesic- 
defined interferometer model. 



estimate the effect due to three nearby gravitational fields: The Earth, the 
Sun, and the Milky Way. 

To use (18. 2p we first need a data point to fix the constant of proportion- 
ality. The data point selected is that with the largest value of R<^. Noting 
that the light travel time differences computed by the reflect routine are in 
units of 2m, where m is the geometric mass of the gravitational field source, 
the data point is 

R^ = 48, = 1, ^ = 2.06 X 10-^ (8.3) 



2m 



From (ESD and (ESD we have 



STre = (2.06 X 10-^)2m^^^^^, (8.4) 

or, since R^, = R/2m, L^ = L/2m, and m = GM/c^ where M is the mass in 
SI units, 

(5rsi = 48^ X (2.06 x 10-^) x (8.5) 

where we have also divided by c to obtain the time difference in seconds, 
rather than metres. 

Using (18. 5p we can estimate the effect due to the Earth, the Sun, and the 
Milky Way. The calculation is summarised in Table 18.11 Compare Table 18.11 
with Table O of Section 17X11 

In Table [Hm the effect due to Milky Way is ~ lO^'^^ s. The smallest time- 
scale currently detectable with gravitational wave detectors is on the order of 
10-2" s. We conclude that the Milky Way cannot be weighed by measuring 

The ordering of the effects (Earth > Sun > Milky Way) is in opposition 
to that of Table 17. 1[ This may be thought of as due to the extra factor of 
L/R in (E2D compared with (El. (For the Milky Way, L/R^ 10-2°, and 
for the Earth, L/i? ~ 10-^.) 
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The estimate for the effect due to the Earth in Table 18.11 cannot be as- 
sumed to be very accurate, because the Schwarzschild radius of the Earth 
is about 9 mm and so, for aim interferometer, L^, ^ 113, which is signifi- 
cantly larger than any value of L^, tested in a numerical experiment in this 
chapter — the relation fl8.2p may not hold in that regime, although we have 
no reason to think it will not. 

8.5 Conclusions 

By simulating the coordinate-defined interferometer of Karim et al. in GR- 
workbench, we were able to validate the theoretical analysis of that interfer- 
ometer, made in [8]. 

By simulating a physically realistic geodesic-defined interferometer, a 
more accurate estimate of the light travel time difference 6Tre on Earth due to 
the Milky Way was obtained, and was found to be too small to be detected. 
It was also found that, in contrast to the case for the coordinate-defined 
interferometer of Karim et ai, the light travel time difference due to the 
gravitational field of the Earth is the most important for an interferometer 
located on Earth, and that due to the gravitational field of the Milky Way is 
the least important of the major gravitational fields in the vicinity of Earth. 

We conclude that the experiment proposed by Karim et al, to weigh the 
galaxy using a small interferometer on Earth, is not feasible, and that their 
conclusion is false because of the approximations implied in their coordinate- 
dependent interferometer model. 



Chapter 9 
Conclusion 



GRworkbench has been successfully and substantially extended to facilitate 
numerical experimentation in General Relativity. 

A functional programming framework has been crucial to the development 
of tools for numerical experimentation within GRworkbench. The functional 
framework is more expressive, permitting important concepts in numerical 
programming and differential geometry to be directly represented in the C++ 
code of GRworkbench. 

New algorithms for key numerical operations have replaced pre-existing 
simpler methods. The numerical engine is now expressed in the paradigm 
of functional programming, enabling algorithms to easily interface with one- 
another. The sophisticated new algorithms are faster and more accurate, and 
an abstraction of the notion of approximate equality enables them to be en- 
coded in a robust and elegant way. Through the C++ template mechanism, 
numerical methods can be encoded such that they can be applied to any sets 
with the required structure defined upon them. 

The differential geometric engine of GRworkbench has been rewritten 
within the functional programming framework. Abstract notions, such as 
points and tangent vectors, are represented by C++ classes. Functions used 
in differential geometry, such as curves in space-time, are now represented 
and manipulated directly as functions. 

Using the new numerical and differential geometric core of GRworkbench, 
tools for numerical experimentation have been developed. Geodesies and 
the parallel transport operation, both implemented in terms of the new ODE 
integration algorithm, represent fundamental physical concepts in General 
Relativity. Methods for determining unique geodesies, defined implicitly in 
terms of boundary conditions, have been developed using the new algorithm 
for function minimisation; these methods enable the construction of photon 
world-lines joining observers to particular events, representing an important 
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physical situation. 

The utihty of numerical experimentation in GRworkbench was demon- 
strated. A traditional analysis of a physical problem in General Relativity, 
involving various simplifying approximations in the mathematical model, was 
investigated and found to yield an inaccurate estimate of the desired phys- 
ical quantity. A more physically motivated model was devised, and an ac- 
curate estimate of the quantity was obtained by simulating the new model 
in GRworkbench. A physically meaningful result was thereby produced by a 
numerical experiment in GRworkbench, where analytic methods had proven 
to be inadequate. 
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Appendix A 

GRworkbench code listings 



This appendix contains code listings from important parts of the rewritten 
numerical engine of GRworkbench, described in Chapter HI and some of the 
tools for numerical experimentation described in Chapter [HI 

Whenever a conflict arises, the algorithms in GRworkbench are generally 
coded with execution speed taking priority over code brevity or simplicity. 
As such, the code in this appendix may appear significantly different to the 
code in Chapters [3l HJ [5l and [61 In many cases, however, the algorithms may 
be more easily read by completely disregarding the symbols const and &l, 
and by interpreting variable declarations of the form 

I const Type variable(expression); 

as the more familiar 

I Type variable = expression; 

A.l Relative difference 

template <typename T> struct relative_difference_implementation 

{ 

static double apply(const J&i a, const J&i b) 

{ 

const double abs_a_abs_b(abs(a) * abs(b)); 
const double abs_a_minus_b(abs(a — b)); 

return abs_a_abs_b <= 1 ? abs_a_minus_b : abs_a_minus_b / sqrt( 
abs_a_abs_b); 

} 

static double apply_squared(const T&: a, const T&i b) 



89 



90 



APPENDIX A. GRWORKBENCH CODE LISTINGS 



{ 

const double abs_a_abs_b(abs(a) * abs(b)); 
const double abs_a_minus_b_squared(square(abs(a — b))); 
return abs_a_abs_b <= 1 ? abs_a_minus_b_squared : 
abs_a_minus_b_squared / abs_a_abs_b; 

} 

}; 

template <typename T> double relative_difference(const a, const T&i 
b) 

{ 

return relative_difference_implementation<T>::apply(a, b); 

} 

template <typename T> double relative_difference_squared(const T&i a, 
const T&i b) 

{ 

return relative_difFerence_implementation<T>::apply_squared(a, b); 

} 

template <typename T> struct relative_difFerence_implementation<nvector 

<T> > 

{ 

static double apply(const nvector<T><Si a, const nvector<T><Si b) 

{ 

return sqrt(apply_squared(a, b)); 

} 

static double apply_squared(const nvector<T>&; a, const nvector<T>&i 
b) 

{ 

if (a.size() != b.size()) 

throw nvector<T>::incompatible(); 

typename nvector<T>::const_iterator i, j; 
double r(0.); 

for (i = a.begin(), j = b.begin(); i != a.end(); ++i, ++j) 
r += relative_difference_squared(*i, *j); 

return r; 

} 

}; 
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template <typename T, size_t N> struct relative_difference_implementation 
<grwb::vector<N, T> > 

{ 

static double apply(const grwb::vector<N, T>&i a, const grwb::vector< 
N, J>ii b) 

{ 

return sqrt(apply_squared(a, b)); 

} 

static double apply_squared(const grwb::vector<N, J>L a, const grwb:: 
vector<N, J>&i b) 

{ 

typename grwb::vector<N, T>::constJterator i, j; 
double r(0.); 

for (i = a.begin(), j = b.begin(); i != a.end(); ++i, ++j) 
r += relative_difference_squared(*i, *j); 

return r; 

} 

}; 

.2 Richardson extrapolation 

template <typename T> class richardson_extrapolation 

{ 

public: 

richardson_extrapolation(const doubled x, const J&i y) 
: limit_(y), 
error_(y) 

{ 

refine(x, y); 

} 

const T&i limit() const 

{ 

return limit.; 

} 

const J&L error() const 
{ 

return error_; 

} 



APPENDIX A. GRWORKBENCH CODE LISTINGS 



void refine(const doubled x, const T&i y) 

{ 

// adapted from Numerical Recipes in C (2nd Edition), p. 731 

data_.resize(data_.size() + 1, make_pair(x, y)); 
error. = limit. = y; 

const size_t n(data_.size()); 
if (n 1) 
return; 

T c(y); 

for (size_t i(l); i < n; ++i) 
{ 

const double x_i(data_[n — i — 1]. first); 

const double delta(l. / (xJ — x)); 

const double fl(x * delta); 

const double f2(x_i * delta); 

const T q(data_[i — 1]. second); 

data_[i — 1]. second = error_; 

const T d2(c - q); 

error. = f 1 * d2; 

c = f 2 * d2; 

limit- += error_; 

} 

data_[n — 1]. second = error_; 

} 

private: 

std::vector<pair<double, T> > data_; 

T limit.; 
T error_; 

}; 

.3 Differentiation 



template <typename T> class derivativeJunctor 
{ 

public: 

derivative_functor(const function<optional<T> (const double&)>& f, 
const doubled scale, const doubled tolerance) 
: f_(f), 



.3. DIFFERENTIATION 



scale_(scale), 
tolerance_(tolerance) 

{ 
} 

optional<T> operator() (const doubled x) const 

{ 

if (!f_(x)) 

return optional<T>(); 

double h(scale_); 

optional<richardson_extrapolation<T> > extrapolator; 

for (size_t i(0); i < max_steps_; ++i) 
{ 

const optional<T> right(f_(x + h)); 
const optional<T> left(f_(x — h)); 
if (left right) 

{ 

const T diff((*right - *left) / (2. * h)); 

if (lextrapolator) 

extrapolator.reset(richardson_extrapolation<T>(h * h, diff)); 
else 

{ 

extrapolator— >refine(h * h, difF); 

if (tolerance. > relative_difFerence(extrapolator— >limit(), 

extrapolator— >limit() + extrapolator— >error())) 
return optional<T>(extrapolator— >limit()); 

} 

} 

h /= step_scale_; 

} 

return optional<T>(); 

} 

private: 

const function<optional<T> (const double(Si)> f_; 
const double scale.; 
const double tolerance.; 
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static const size_t max_steps_ = 13; 
static const double step_scale_ = 1.7; 

}; 

template <typename T> function<optional<T> (const double&)> 

derivative(const function<optional<T> (const double<Si)>&i f, const 
double&i scale = 1., const double&i tolerance = 
default_approx_equaLtolerance) 

{ 

return derivative_functor<T>(f, scale, tolerance); 

} 

A. 3.1 Gradient 

namespace gradient_detail 

{ 

template <typename T> class single_coordinate_function 

{ 

public: 

single_coordinate_function(const function<optional<T> (const nvector 
<double>&)><£i f, const nvector<double>(S^ x, const s\zeJ.&i i) 
: f_(f), 
x-(x), 
i.(i) 

{ 
} 

optional<T> operator() (const doubled delta_xJ) 

{ 

nvector<double> _x(x_); 
_x[i_] += delta_x:_i; 
return f-(-x); 

} 

private: 

const function<optional<T> (const nvector<double>&)>.Si f_; 
const nvector<double>&; x_; 
const size_t<^ i_; 

}; 

} 



template <typename T> class gradient_functor 
{ 
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public: 

gradient_functor(function<optional<T> (const nvector<double>&)>f) 
: f.(f) 

{ 
} 

optional<nvector<T> > operator() (const nvector<double>.Si x) 

{ 

const optional<T> default_value(f_(x)); 

if (!default_value) 

return optional<nvector<T> >(); 

nvector<T> result(x.size(), unchanging(*default_value)); 

for (size_t i = 0; i != x.size(); 

{ 

optional<T> d(derivative<T>(gradient_detail:: 

single_coordinate_function<T>(f_, x, i))(0.)); 
if (!d) 

return optional<nvector<T> >(); 
result[i] = *d; 

} 

return optional<nvector<T> >(result); 

} 

private: 

const function<optional<T> (const nvector<double>&i)> f_; 

}; 

template<typename T> function<optional<nvector<T> >(const nvector< 
double>&i)> gradient(const function<optional<T>(const nvector< 
double>&)>& f) 

{ 

return gradient_functor<T>(f); 

} 

A. 4 Bulirsch-Stoer method 

template <class T, template <class> class U> class bulirsch_stoer 

// adapted from Numerical Recipes in C (2nd Edition), p. 728 

{ 
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public: 

bulirsch_stoer(const function<optional<T> (const doubled, const Tii) 
>&i f, const double&i x_0, const y_0, const double&i 
default-Stepsize = 1., const size_t&; maximum_steps = 100, const 
doubled relative_error = default_approx_equaLtolerance) 
: f.(f), 

maximum_steps_(maximum_steps), 
relative_error_(relative_error), 
defa u lt_h_(defa u lt_stepsize) , 
x_(x_0), 

y-(y-o) 



{ 



const double safe_relative_error(relative_error * safel_); 

typename std::map<double, vector<bulirsch_stoer_parameters<U>:: 
k_total, vector<bulirsch_stoer_parameters<U>::k_total, double> > 
> : :const_iterator i (a I pha_cache_() .f i nd (safe_relative_error) ) ; 

if (i != alpha_cache_().end()) 
alpha. = i— >second; 

else 

{ 

for (size_t i = 1; i < bulirsch_stoer_parameters<U>::k_total; ++i) 
for (size_t j = 0; j < i; ++j) 

alpha_p][i] = pow(safe_relative_error, (a_()[j] — a_()[i]) / ((a. 

()[i]-a_()[0] + l.)*(2*j + 3))); 
alpha_cache_()[safe_relative_error] — alpha.; 

} 

for (optimaLk_ = 1; optimal_k_ < bulirsch_stoer_parameters<U>;:k_total 
— 1; ++optimaLk_) 

if (a_()[optimaLk_ + 1] > a_()[optimaLk_] * alpha_[optimaLk_ — 1][ 
optimal_k_]) 
break; 
max_k_ = optimaLk_; 



} 



const double<^ x() const 

{ 

return x_; 

} 



const y() const 

{ 
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return y_; 

} 

bool step(const double&i to_x) 
{ 

double h(default_h_); 
if (to_x < x_) 
h *= -1; 

for (size_t i = 0; i < maximum_steps_; ++i) 

{ 

bool reduced_step_size(false); 

bool success(false); 

size_t k(0), km(0); 

double stepsize_reduction_factor(0.); 

double err[bulirsch_stoer_parameters<U>::k_total]; 

U<T> stepper(f_, x_, y_); 

if (to_x == x_) 
return true; 

if ((to_x - x_) * (to_x - x_ - h) < 0.) 
h — to_x — x_; 

while (true) 

{ 

optional<richardson_extrapolation<T> > extrapolator; 

for (k = 0; k < max_k_; ++k) 
{ 

optional<T> y_est(stepper.step(h, bulirsch_stoer_parameters<U 

>::k_values[k])); 
if (!y-est) 

return false; 

const double little_h_squared(square(h / 

bulirsch_stoer_parameters<U>::k_values[k])); 

if (lextrapolator) 

extra polator.reset(richardson_extrapolation<T>( 
little_h_squared, *y_est)); 

else 

{ 

extrapolator— >refine(little_h_squared, *y_est); 
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y_ = extrapolator— >limit(); 

const double error(relative_difference(y_, y_ + extrapolator 

— >error()) / relative_error_); 
km = k — 1; 

err[km] = pow(error / safel_, 1. / (2 * km + 3)); 

if (k >= optimaLk. - 1 || i == 0) 
{ 

if (error < 1.) 

{ 

success = true; 
break; 

} 

if (k max_k_ || k optimaLk_ + 1) 

{ 

stepsize_reduction_factor = safe2_ / err[km]; 
break; 

} 

if (k == optimaLk_ alpha_[optimal_k_ — 1][ 
optimaLk_] < err[km]) 

{ 

stepsize_reduction_factor — 1- / err[km]; 
break; 

} 

if (optimaLk_ —— max_k_ alpha_[km][max_k_ — 1] < 
err[km]) 

{ 

stepsize_reduction_factor = alpha_[km][max_k_ — 1] 

* safe2_ / err[km]; 
break; 

} 

if (alpha_[km][optimal_k_] < err[km]) 
{ 

stepsize_reduction_factor — alpha_[km][optimaLk_ — 

1] / err[km]; 
break; 

} 

} 

} 

} 



if (success) 
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break; 

if (stepsize_reduction_factor > min_stepsize_reduction_) 
stepsize_reduction_factor = min_stepsize_reduction_; 
if (stepsize_reduction_factor < max_stepsize_reduction_) 
stepsize_reduction_factor = max_stepsize_reduction_; 
h *= stepsize_reduction_factor; 
reduced_step_size = true; 

} 

x_ += h; 

double work_min(l.e300); 

double scale_factor(0.); 

for (size_t j = 0; j <= km; ++j) 

{ 

const double s(err[j] < max_stepsize_increase_ ? 

max_stepsizeJncrease_ : errp]); 
const double work(s * a_()[j + 1]); 
if (work < work_min) 

{ 

scale_factor — s; 
work_min = work; 
optimaLk_ = j + 1; 

} 

} 

if (optimaLk_ >= k optimaLk_ != max_k_ !reduced_step_size) 

{ 

double s(scale_factor / alpha_[optimaLk_ — l][optimaLk_]); 
if (s < max_stepsizeJncrease_) 
s — max_stepsizeJncrease_; 

if (a_()[optimal_k_ + 1] * s <= work_min) 
{ 

sea I e_f actor = s; 
++optimaLk_; 

} 

} 

h /= scale_factor; 

} 

cout << " Bulirsch— Stoer:^Too^many^teps^required." << endl; 
return false; 
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} 

private: 

const static double safel_ = 0.25; 

const static double safe2_ = 0.7; 

const static double max_stepsize_reduction_ = l.e— 5; 

const static double min_stepsize_reduction_ = 0.7; 

const static double max_stepsize_increase_ = 0.1; 

const function<optional<T> (const doubled, const Tii)> f_; 
const size_t maximum_steps_; 
const double relative_error_; 
const double default-h_; 

vector< bu 1 i rsch_stoer_para meters< U > : : k_tota I , vector< 

bulirsch_stoer_parameters<U>::k_total, double> > alpha.; 

static std::map<double, vector<bulirsch_stoer_parameters<U>::k_total, 
vector<bulirsch_stoer_parameters<U>::k_total, double> > 
alpha_cache_() 

{ 

static std::map<double, vector<bulirsch_stoer_parameters<U>::k_total, 

vector<bulirsch_stoer_parameters<U>::k_total, double> > > _; 
return _; 

}; 

static vector<bulirsch_stoer_parameters<U>::k_total + 1, double>& a_() 

{ 

static optional<vector<bulirsch_stoer_parameters<U>::k_total + 1, 

double> > _; 
if (L) 

{ 

_.reset(vector<bulirsch_stoer_parameters<U>::k_total + 1, double>()); 

(*_)[0] = bulirsch_stoer_parameters<U>::k_values[0] + 1; 

for (size_t i = 0; i < bulirsch_stoer_parameters<U>::k_total; ++i) 

(*_)[i + 1] = (*-)[i] + bulirsch_stoer_parameters<U>::k_values[i + 

1]; 

} 

return *_; 

} 

double x_; 

Ty-; 
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size_t optimaLk_; 
size_t max_k_; 

}; 

A. 4.1 Modified midpoint method 

template <class T> class mociified_midpoint_stepper 
{ 

public: 

modified_midpoint_stepper(const function<optional<T> (const doubled, 
const T&i)>&i f, const double<S^ x_0, const J&i y_0) 
: f.(f), 
x_0_(x_0), 

y-O-(y-O), 

f_y_0_(f(x_0, y_0)) 

{ 
} 

optional<T> step(const doubled totaLh, const size_t&i steps) const 
{ 

// adapted from Numerical Recipes in C (2nd Edition), p. 724 

optional<T> ret; 

if (!f-y-0_) 
return ret; 

const double h(totaLh / double(steps)); 
const double two_h(2. * h); 

T ym(y_0_), yn(y_0_ + h * *f_y_0_); 
double x(x_0_ + h); 

optional<T> dydx(f_(x, yn)); 
if (!dydx) 
return ret; 

for (size_t i(l); i < steps; 

{ 

T y_next(ym + two_h * *dydx); 
ym = yn; 

dydx = f_(x += h, yn = y_next); 
if (!dydx) 
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return ret; 

} 

ret.reset(0.5 * (ym + yn + h * *dyclx)); 
return ret; 

} 

private: 

const function<optional<T> (const doubled, const Tii)> f_; 

const double x_0-; 

const T y_0_; 

const optional<T> f_y_0_; 

}; 

template <> class bulirsch_stoer_parameters<modified_midpoint_stepper> 
{ 

public: 

const static size_t k_total — 10; 

const static size_t k_values[k_total + 1]; 

}; 

template <> const size_t bulirsch_stoer_parameters< 

modified_midpoint_stepper>::k_values[] = {2, 4, 6, 8, 10, 12, 14, 16, 18, 
20, 22}; 

A. 5 Powell's method 



template <typename T, typename U> class powelLmini miser 

{ 

public: 

powelLminimiser(const function<optional<T> (const U&i)>&; f) 
: f_(f) 

{ 
} 

optional<pair<U, T> > operator() (const U& x, const nvector<U>.Si 
basis, const doubled tolerance = default_approx_equaLtolerance) 
const 

{ 

U minimum(x); 

optional<T> op(f_(minimum)); 
if (lop) 

return optional<pair<U, T> >(); 
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T f_min(*op); 
nvector<U> basis_( basis); 

for (size_t i(0); i < max_steps; ++i) 

{ 

const U prev_min(minimum); 
const T prev_f_min(f_min); 
T largest_decrease(zero(f_min)); 
size_t largest_decrease_index(0); 

for (size_t j(0); j < basis.size(); ++j) 
{ 

optional<pair<double, T> > line_minimum(brent_minimiser( 
linear_subspace(f_, minimum, basis_p]))(0., 0., tolerance)); 

if (!line_minimum) 

return optional<pair<U, T> >(); 

if (f_min — line_minimum— >second > largest_decrease) 
{ 

largest-decrease = f_min — line_minimum— >second; 
largest_decrease_index = j; 

} 

if (zero(line_minimum— >first) != Iine_minimum— >first f_min > 
line_minimum— >second) 

{ 

basis_[j] *= line_minimum— >first; 

minimum += basis_[j]; 

f_min = line_minimum— >second; 

} 

} 

if (approx_equal(prev_f_min, f_min, tolerance)) 

return optional<pair<U, T> >(make_pair(minimum, f_min)); 

const U new_direction(minimum — prev_min); 
const U extrapolated_min(minimum + new_direction); 
op — f_(extrapolated_min); 
if (lop) 

return optional<pair<U, T> >(); 
const T f_extrapolated_min(*op); 

if (f_extrapolated_min < prev_f_min) 
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{ 

if (2. * (prev_f_min — 2. * f_min + f_extrapolated_min) * square( 
prev_f_min — f_min — f_extrapolated_min) <= largest_decrease * 
square(prev_f_min — f_extrapolated_min)) 

{ 

optional<pair<double, T> > line_minimum(brent_minimiser( 
linear_subspace(f_, minimum, new_direction))(0., 0., tolerance 

)); 

if (!line_minimum) 

return optional<pair<U, T> >(); 

if (zero(line_minimum— >first) != Iine_minimum— >first f_min 
> line_minimum— >second) 

{ 

basis_[largest_decrease_index] = basis_[basis.size() — 1]; 
minimum += (basis_[basis.size() — 1] = new_direction * 

line_minimum— >first); 
f_min — line_minimum— >second; 

} 

} 

} 

} 

return optionai<pair<U, T> >(); 

} 

optional<pair<U, T> > operator()(const U& x, const doubled tolerance 
= default_approx_equaLtolerance) const 

{ 

return operator()(x, default_basis(x), tolerance); 

} 

private: 

const function<optional<T> (const U&)> f_; 

const static size_t max_steps = 100; 
const static double auto_scale = l.e— 2; 

class linear_subspace_functor 

{ 

public: 

linear_subspace_functor(const function<optional<T> (const U&)><Si f, 
const Uii origin, const \Jii direction) 
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: f_(f), 
origin_(origin), 
direction .(direction) 

{ 
} 

optional<T> operator() (const doubled t) const 

{ 

return f_(origin_ + t * direction.); 

} 

private: 

const function<optional<T> (const U8i)>&i f_; 
const U& origin.; 
const U&; direction.; 

}; 

function<optional<T> (const double<Si)> linear.subspace(const function< 
optional<T> (const U&)>& f, const origin, const U&i direction) 
const 

{ 

return linear_subspace_functor(f, origin, direction); 

} 

nvector<U> default_basis(const U&i x) const 

{ 

nvector<U> r(unity(nvector<U>(x.size(), unchanging(x)))); 
for (size_t i(0); i < x.size(); ++i) 

{ 

const double scale(auto_scale * abs(x[i])); 
if (scale > 0) 
r[i] *= scale; 

} 

return r; 

} 

}; 

A. 5.1 Brent minimiser 

template <typename T, typename U> class bracketer 
{ 

public: 

bracketer(const function<optionai<T> (const U<S^)>&; f) 
:_(f) 
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{ 
} 

optional<vector<3, pair<U, T> > > operator() (const \J&i x, const \J&i 
step_size — U()) const 

{ 

optional<vector<3, pair<U, T> > > result; 

optional<T> op(_(x)); 
if (!op) 
return result; 

vector<3, pair<U, T> > r(unchanging(make_pair(x, *op))); 
U step = step_size —— U() ? (x zero(x) ? unity(x) : auto_scale * abs( 
x)) : step_size; 

r[l]. first += step; 
op = _(r[l].first); 
if (lop) 

return result; 
r[l]. second = *op; 

if (r[l]. second > r[0]. second) 

{ 

swap(r[0], r[l]); 
step *= — 1; 

} 

for (size_t i(0); i < max_steps; ++i) 
{ 

r[2]. first = r[l]. first + step; 
op = _(r[2].first); 
if (lop) 

return result; 
r[2]. second = *op; 

if (r[2]. second >= r[l]. second) 
{ 

if (step < 0) 

swap(r[0], r[2]); 
result. reset(r); 
return result; 

} 
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r[0] = r[l]; 
r[l] = r[2]; 
step ratio; 

} 

return result; 

} 

private: 

const function<optional<T> (const \Jii)> _; 

const static size_t max_steps = 100; 

const static double ratio = 1.6; 

const static double auto_scale = l.e— 2; 

}; 

template <typename T, typename U> class brent_minimiser_functor 
// adapted from Numeric Recipes in C (2nd Edition), p. 404 

{ 

public: 

brent_minimiser_functor(const function<optional<T> (const U&)>&i f) 
: f.(f) 

, bracketer_(f) 

{ 
} 

optional<pair<U, T> > operator() (const U& x, const U&: scale = U(), 
const double&i tolerance = default_approx_equaLtolerance) const 

{ 

const optional<vector<3, pair<U, T> > > bracket(bracketer_(x, scale)); 
if (! bracket) 
return optional<pair<U, T> >(); 

U left((*bracket)[0]. first), best((*bracket)[l].first), right((*bracket)[2].first 

); 

U third_best(best), second_best(best); 
T f_best((*bracket)[l]. second); 

T f_trial(f_best), f_third_best(f_best), f_second_best(f_best); 
U d(0.), prev_d(0.); 

const double two_tolerance(2. * tolerance); 
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for (size_t i(0); i < max_steps; ++i) 
{ 

const U mid(0.5 * (left + right)); 

if (approx_equal(left, right, two_tolerance)) 
return optional<pair<U, T> >(make_pair(best, f_best)); 

if (abs(prev_d) > tolerance * abs(best)) 
{ 

const U r((best — second-best) * (f_best — f_third_best)); 

U q((best — third_best) * (Lbest — f_second_best)); 
U p((best — third_best) * q — (best — second-best) * r); 
q = 2. * (q - r); 
if (q > 0.) 

p = -p; 

q = abs(q); 

U prev_prev_d(prev_d); 
prev_d = d; 

if (abs(p) >= abs(0.5 * q * prev_prev_d) || p <= q * (left — best) || 
p >= q * (right — best)) 

d = cgold * (prev_d — (best >= mid ? left — best : right — best 

)); 

else 

d = P / q; 

} 

else 

d = cgold * (prev_d = best >= mid ? left — best : right — best); 

const U trial(best + d); 

const optional<T> op(f_(trial)); 

if (lop) 

return optional<pair<U, T> >(); 
f_trial = *op; 

if (f_trial <= f_best) 
{ 

if (trial >= best) 
left = best; 
else 

right = best; 

third_best = second-best; second-best = best; best = trial; 
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f_third_best — f_second_best; f_second_best = f_best; Lbest = Ltrial 

} 

else 

{ 

if (trial < best) 
left = trial; 

else 

right = trial; 

if (Ltrial <= f_second_best || second-best —— best) 
{ 

third_best = second_best; second_best = trial; 
f_third_best = f_second_best; f_second_best = Ltrial; 

} 

else if (Ltrial <= Lthird_best || third_best best || third_best —- 
second-best) 

{ 

third_best = trial; 
Lthird_best = Ltrial; 

} 

} 

} 

return optional<pair<U, T> >(); 

} 

private: 

const function<optional<T> (const \Jii)> L; 
const bracketer<T, U> bracketer_; 

const static size_t max_steps = 100; 
const static double cgold — 0.3819660; 

}; 

template <class T, class U> brent_minimiser_functor<T, U> 

brent_minimiser(const function<optional<T> (const U<Si)>& f) 

{ 

return brent_minimiser_functor<T, U>(f); 

} 

A. 6 Geodesic 



class geodesic 
{ 
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public: 

geodesic(const tangent_vector& t) 
: atlas_(t.context().context()), 
least_upper_bou nd_(positi veJ nf i n ity) , 
greatestJower_bound_(negative_infinity), 
cache_(new std::map<double, tangent_vector>) 

{ 

cache.— >insert(cache_value_type(0, t)); 

} 

optional<point> operator() (const doubled t) const 

{ 

return operator()(t, 0.). second; 

} 

pair<double, optional<point> > operator()(const doubled t, const 
doubled epsilon) const 

{ 

typedef pair<double, optional<point> > return_type; 

if (t >= least_upper_bound_) 

return return_type(least_upper_bound_, optional<point>()); 
else if (t <= greatest_lower_bound_) 

return return_type(greatestJower_bound_, optional<point>()); 

cache_iterator_type initiaLdata(get_initiaLdata(t)); 

if (abs(initiaLdata— >first — t) <= epsilon) 
return return_type(initial_data— >first, optional<point>(initiaLdata— > 
second. context())); 

optional<cache_iterator_type> result(advance(initiaLdata— >second, 
initiaLdata— >first, t)); 

if (result) 

return return_type(t, optional<point>((*result)— >second.context())); 

if (t > 0) 

least_upper_bound_ — t; 
else 

greatestJower_bound_ = t; 
return return_type(t, optional<point>()); 
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} 

optional<tangent_vector> tangent(const doubled t) const 

{ 

if (!operator()(t)) 
return optional<tangent_vector>(); 

return optional<tangent_vector>(cache_— >find(t)— >second); 

} 

private: 
const weak_ptr<atlas> atlas_; 

mutable double least-upper.bound.; 
mutable double greatest_lower_bound_; 

shared_ptr<std::map<double, tangent_vector> > cache.; 

typedef std::map<double, tangent_vector>::const_iterator 

cache_iterator_type; 
typedef std::map<double, tangent_vector>::value_type cache_value_type; 

class geodesic-callback 

{ 

public: 

geodesic_callback(const shared_ptr<atlas::chart>(S^ c) 
:.(c) 

{ 
} 

optional<vector<2, nvector<double> > > operator() (const doubled t, 
const vector<2, nvector<double> >&; y) const 

{ 

const int dim(y[0].size()); 

vector<2, nvector<double> > ret(unchanging(y[l])); 

optional<ntensor_components<3>::type> con(connection(*_)(y[0])); 
if (Icon) 

return optional<vector<2, nvector<double> > >(); 
ret[l] *= 0.; 

for (int a = 0; a < dim; ++a) 
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for (int b = 0; b < dim; ++b) 
for (int c = 0; c < dim; ++c) 

ret[l][a] -= (*con)[a][b][c] * y[l][b] * y[l][c]; 

return optional<vector<2, nvector<double> > >(ret); 

} 

private: 

const shared_ptr<atlas::chart> _; 

}; 

cache_iterator_type get_initial_data(const double&i t) const 

{ 

cache_iterator_type after(cache_— >lower_bound(t)); 

if (after == cache.— >begin()) 
return after; 

if (after == cache_— >end()) 
return — cache_iterator_type(after); 

cache_iterator_type before(after); 

before; 

return abs(before— >first — t) < abs(after— >first — t) ? before : after; 

} 

optional<cache_iterator_type> advance(const tangent_vector(S^ tangent, 
const double<S^ from_t, const double&i to_t, size_t recursion = 1) 
const 

{ 

const pointiSi origin — tangent. context(); 

if (recursion > max_recursion_) 
return optional<cache_iterator_type>(); 

//cout « " Geodesic: " « from.t «"->"« to.t « end I; 

shared_ptr<atlas> a(atlas_); 

for (set<shared_ptr<atlas::chart> >::const_iterator i = a— >charts.begin 
(); i != a->charts.end(); ++i) 
if (origin[*i] tangent[*i]) 

{ 
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optional<cache_iterator_type> result(advance_on_chart(*i, *origin[*i], 

*tangent[*i], from_t, to_t)); 
if (result) 

return result; 

} 

optional<cache_iterator_type> halfway(advance(tangent, from_t, (from_t 

+ to_t) / 2, recursion + 1)); 
return halfway ? advance((*halfway)— >second, (from_t + to_t) / 2, to_t, 
recursion + 1) : optional<cache_iterator_type>(); 

} 

optional<cache_iterator_type> advance_on_chart(const shared_ptr<atlas:: 
chart>& c, const nvector<double>& x, const nvector<double><£i dx, 
const doubled from.t, const doubled to.t) const 

{ 

bulirsch_stoer<vector<2, nvector<double> >, modified_midpoint_stepper 
> solver((geodesic_callback(c)), from_t, make_vector(x, dx)); 

if (!solver.step(to_t)) 
return optional<cache_iterator_type>(); 

point dest(atlas_.lock(), c, solver. y()[0]); 
tangent_vector tv(dest, c, solver. y()[l]); 

return optional<cache_iterator_type>(cache_— >insert(cache_value_type( 
to_t, tv)). first); 

} 

const static size_t max_recursion_ — 7; 

const static double positiveJnfinity = leSOO; 
const static double negativeJnfinity = — leSOO; 

}; 

A. 7 Generalised spherical polar coordinate trans- 
formation 

template <typename T> inline nvector<T> to_polar(const nvector<T>& 
x) 

{ 

nvector<T> r(x); 
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T cosine(abs(x)); 

const T zero_t(zero(cosine)); 

r[0] = cosine; 

for (size_t i(l); i < x.size() - 1; ++i) 
if (cosine != zero_t) 

{ 

r[i] = asin(x[i — 1] / cosine); 
cosine cos(r[i]); 

} 

else 

r[i] = zero_t; 

r[x.size() — 1] = atan2(x[x.size() — 2], x[x.size() — 1]); 
return r; 

} 

template <typename T> inline nvector<T> from_polar(const nvector<T 

>&L X) 

{ 

nvector<T> r(x); 
T cosine(x[0]); 

for (size_t 1(0); i < x.size() - 1; ++i) 
{ 

r[i] = cosine * sin(x[i + 1]); 
cosine *= cos(x[i + 1]); 

} 

r[x.size() — 1] = cosine; 
return r; 

} 

template <typename T> inline nvector<T> from_polar_with_radius(const 
nvector<T>(^ x, const J&i radius) 

{ 

nvector<T> x2(x.size() + 1, unchanging(radius)); 
for (size_t i(0); i < x.size(); ++i) 

x2[i + 1] = x[i]; 
return from_polar(x2); 

} 
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template <typename T> inline nvector<T> to_polar_without_radius(const 
nvector<T>& x) 

{ 

const nvector<T> polar(to_polar(x)); 

return nvector<T>(polar.size() — 1, polar. begin() + 1); 

} 

A. 8 Connecting geodesic 

namespace connecting_geodesic_cletail 
{ 

class geodesic-shooter 

{ 

public: 

geodesic_shooter(const point&c a, const points b, const shared_ptr<atlas 
::chart><Si c) 
: a_(a), 
b-(b), 
chart_(c) 

{ 
} 

optional<double> operator()(const nvector<double>&; v) const 

{ 

geodesic geo(tangent_vector(a_, chart., from_polar_with_radius(v, 1.))); 
optional<pair<double, double> > r(min_euclidean_separation(geo, b_) 

); 

return r ? optional<double>(r— >second) : optional<double>(); 

} 

private: 
const point&i a_; 
const point&: b_; 

const shared_ptr<atlas::chart>&i chart.; 

}; 

} 

inline optional<geodesic> connecting_geodesic(const point<S^ a, const point<S^ 
b) 

{ 

const shared_ptr<atlas::chart>&i c(a.valid_chart()); 
connecting_geodesic_detaih:geodesic_shooter shooter(a, b, c); 
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function<optional<double> (const nvector<double>&;)> 

shooter_f u nction (shooter) ; 
powelLminimiser<double, nvector<double> > pm(shooter_function); 

const nvector<double> va(*a[c]), vb(*b[c]); 

optional<pair<nvector<double>, double> > r(pm(to_polar_without_radius 
(vb - va))); 

if (!r) 

return optional<geodesic>(); 

const nvector<double> v(from_polar_with_radius(r— >first, 1.)); 
const double scale(min_euclidean_separation(geodesic(tangent_vector(a, c, v 
)), b)->first); 

return optional<geodesic>(geodesic(tangent_vector(a, c, v * scale))); 

} 

A. 9 Connecting null geodesic 

namespace connecting_null_geodesic_detail 

{ 

class null_geodesic_shooter 

{ 

public: 

nulLgeodesic_shooter(const function<optional<point> (const double&i) 
>&i curve, const point<S^ a, const shared_ptr<atlas::chart>&; c, 
const nvector<nvector<double> ><S^ tangent_basis, const doubled 
nulLguess, const doubled guess) 
: guess_(make_vector(nulLguess, guess)), 

curve_(curve), 

a-(a). 
chart_(c), 

basis(tangent_basis) 

{ 
} 

optional<double> operator() (const nvector<double>iS^ v) const 

{ 

const nvector<double> v2(from_polar_with_radius(v, 1.)); 
nvector<double> v3(v2.size() + 1, unchanging(l.)); 
for (size.t i(0); i < v2.size(); ++i) 
v3[i + 1] = v2[i]; 
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geodesic geo(tangent_vector(a_, chart., basis * v3)); 
optional<pair<nvector<double>, double> > r( 
min_euclidean_separation(geo, curve., guess_)); 

if (!r) 

return optional<double>(); 

guess_ = r— >first; 

return optional<double>(r— >second); 

} 

const nvector<double><Si last_guess() const 
{ 

return guess_; 

} 

private: 

mutable nvector<double> guess_; 

const function<optional<point> (const double&;)>&i curve.; 
const point&i a_; 

const shared_ptr<atlas::chart><Si chart.; 
const nvector<nvector<double> basis; 

}; 

} 

inline optional<pair<double, geodesio > connecting_null_geodesic(const 
function<optional<point> (const double&;)>(Si curve, const point<Si a, 
const double&i guess) 

{ 

const static nvector<nvector<double> > polar_basis(make_vector( 
make_vector(0.01, 0.), make_vector(0., 0.01))); 

const shared_ptr<atlas::chart><Si c(a.vaiid_chart()); 

const nvector<nvector<double> > basis(orthonormaLtangent_basis(a, c)); 

const point b(*curve(guess)); 

const nvector<double> va(*a[c]), vb(*b[c]); 

const nvector<double> spacelike(vb.size() — 1, (inverse(basis) * (vb — va) 
).begin() + 1); 

const nvector<double> spacelike_polar(to_polar_without_radius(spacelike)); 
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const connecting_nulLgeodesic_detail::null_geodesic_shooter shooter(curve, a 

, c, basis, abs(spacelike), guess); 
const function<optional<double> (const nvector<double>&)> 

shooterJunction (shooter); 
const powelLminimiser<double, nvector<double> > pm(shooter_function) 

optional<pair<nvector<double>, double> > r(pm(spacelike_polar, 

polar_basis)); 
if(!r) 

return optional<pair<double, geodesio >(); 

const nvector<double> vr(from_polar_with_radius(r— >first, 1.)); 
nvector<double> vr2(vr.size() + 1, unchanging(l.)); 
for (size_t i(0); i < vr.size(); ++i) 

vr2[i + 1] = vr[i]; 
const nvector<double> solution(basis * vr2); 

const nvector<double> scales(min_euclidean_separation(geodesic( 

tangent_vector(a, c, solution)), curve, sliooter.last_guess())— >first); 

return optional<pair<double, geodesio >(make_pair(scales[l], geodesic( 
tangent_vector(a, c, solution * scales[0])))); 

} 

inline optional<pair<double, geodesio > connecting_nulLgeodesic(const 
point&: a, const function<optional<point> (const double&)>& curve, 
const double&i guess = 0.) 

{ 

return connecting_nulLgeodesic(curve, a, guess); 

} 



Appendix B 

Numerical experiment code 
listing 

This appendix lists the code, written by the author, for the numerical ex- 
periment described in Chapter [HI The code fragments listed in that chapter 
were adapted from portions of this code; they were simplified for clarity, and 
any code not directly relevant to the discussion was removed. 

The comments at the beginning of Appendix |X] also apply to the code 
listed here. 

template <class T> class geodesicJnterferometer ; public 
numericaLexperiment<T> 

{ 

public: 

explicit geodesic_interferometer(const shared_ptr<T>&; _); 
private: 

void reflect(const doubled r, const doubled interferometer_speed, const 
double&i arm_length); 

}; 

template <class T> class karimJnterferometer : public 
numerical_experiment<T> 

{ 

public: 

explicit karim_interferometer(const shared_ptr<T>& _); 
private: 

void reflect(const double&i r, const doubled interferometer_speed, const 
doubled arm_length); 

}; 
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template <class T> inline karim_interferometer<T>::karim_interferometer( 
const shared_ptr<T><?^ _) 
: numerical_experiment<T>(_) 

{ 

add_distortion(" SphericaLto^Orthonormal" ); 

add_distortion(" Linear" ); 

nvector<nvector<double> >&; lin(*dynamic_pointer_cast<linear_distortion 

>(distortions.back())— >matrix); 
lin[0][0] = 0.; 
Iin[0][3] = 1.; 

reflect(10., 0.2, 3.); 

} 

template <class T> inline geodesic_interferometer<T>:: 
geodesic_interferometer(const shared_ptr<T>&i _) 
: numericaLexperiment<T>(_) 

{ 

add_distortion("Spherical^to^Orthonormal"); 
add_distortion(" Linear" ); 

nvector<nvector<double> lin(*dynamic_pointer_cast<linear_distortion 

> (distortions. back())—>matrix); 
lin[0][0] = 0.; 
Iin[0][3] = 1.; 

reflect(4., 0.2, 2.); 

} 

template <class T> inline void geodesic_interferometer<T>::reflect(const 
doublec^ r, const doubled interferometer_speed, const doubled 

armJength) 

{ 

const double two_arm_length(2. * armJength); 

const shared_ptr<atlas::chart>&; c(any_chart()); 

const nvector<double> coordinate_direction(make_vector(l., 0., 0., 
interferometer_speed / r)); 

const point origin(this— >atlas(), c, make_vector(0., r, halLpi, halLpi)); 
const tangent_vector tangent(normalise(tangent_vector(origin, c, 
coordinate_direction))); 
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const cached_worldline beam_splitter(coordinate_line(tangent, c)); 

cout << endl << " Beam— splitter^origin^=J' << *origin[c] << ".Jangent 

^—J « *tangent[c] << "." << endl; 
cout << " lnterferometer^peed:J' << interferometer_speed << endl; 
cout << " lnterferometer^armJength:J' << armJength << "." << endl; 
output << r << endl; 
output << interferometer_speed << endl; 
output << armJength << endl; 

const tangent_vector radiaLmirror_direction(orthonormalise(tangent_vector( 

origin, c, make_vector(0., 1., 0., 0.)), tangent)); 
const tangent_vector theta_mirror_direction(orthonormalise(orthonormalise( 

tangent_vector(origin, c, make_vector(0., 0., 1., 0.)), tangent), 

radiaLmirror_direction)); 
const tangent_vector phi_mirror_direction(orthonormalise(orthonormalise( 

orthonormalise(tangent_vector(origin, c, make_vector(0., 0., 0., 1.)), 

tangent), radiaLmirror_direction), theta_mirror_direction)); 
const vector<5, tangent_vector> directions(make_vector( 

radiaLmirror_direction, — radiaLmirror_direction, theta_mirror_direction, 

phi_mirror_direction, — phi_mirror_direction)); 

double max_affine_length(0.); 

for (const tangent_vector* i(directions.begin()); i != directions. end(); ++i) 
{ 

const point mirror_origin(*geodesic(*i)(arm_length)); 

const tangent_vector mirror_tangent(normalise(tangent_vector( 

mirror_origin, c, coordinate_direction))); 
const cached_worldline mirror(coordinateJine(mirror_tangent, c)); 

const geodesic outward_ray(connecting_null_geodesic(origin, mirror, 

armJength)— >second); 
const point reflection(*outward_ray(l.)); 

const pair<double, geodesic> con(*connecting_nulLgeodesic(reflection, 

beam_splitter, two_armJength)); 
const doubled affineJength(con.first); 
const geodesiciSi inward_ray(con. second); 

plot(mirror, 0., affineJength); 
plot(outward_ray, 0., 1.); 
plot(inward_ray, 0., 1.); 
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cout << "^^Mirror^irection^=J' << *(*i)[c] << " , Jime^experienced^ 
by^beam^plitter^=J' << affineJength << "." << endl; 

cout << "^^^^Outward jay^tangent^=J' << *(*outward_ray.tangent(0.)) 
[c] << endl; 

cout << "^^^Jnward^rayJangent^=J' << *(*inward_ray.tangent(0.))[c] 

<< endl; 
output << affineJength << endl; 

if (affineJength > max_affineJength) 
max_affineJength = affineJength; 

} 

plot(beam_splitter, 0., max_affineJength); 

} 

template <class T> inline void karimJnterferometer<T>::reflect(const 
doubled r, const doubled interferometer_speed, const doubled 

armJength) 

{ 

const double two_arm_length(2. * armJength); 

const shared_ptr<atlas::chart>& c(any_chart()); 

const nvector<double> coordinate_direction(make_vector(l., 0., 0., 
interferometer_speed / r)); 

const point origin(this— >atlas(), c, make_vector(0., r, half_pi, half_pi)); 
const tangent_vector tangent(normalise(tangent_vector(origin, c, 

coordinate_direction))); 
const cached_worldline beam_splitter(coordinateJine(tangent, c)); 

cout << endl << " Beam— splitterjjrigin: J' << *origin[c] << ".^tangent:^ 

<< *tangent[c] << "." << endl; 
cout << " Interferometer^peed: J' << interferometer_speed << endl; 
cout << " Interferometer^armJength: J' << armJength << "." << endl; 
output << r << endl; 
output << interferometer_speed << endl; 
output << armJength << endl; 

const tangent_vector radial_mirror_direction(origin, c, make_vector(0., 1., 0., 

0.)); 

const tangent_vector theta_mirror_direction(origin, c, make_vector(0., 0., 1. 
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/ r, 0.)); 

const tangent_vector phi_mirror_direction(origin, c, make_vector(0., 0., 0., 1. 

/O); 

const vector<5, tangent_vector> directions(make_vector( 

radiaLmirror_direction, — radiaLmirror_direction, theta_mirror_direction, 
phi_mirror_direction, — phi_mirror_direction)); 

double max_affineJength(0.); 

for (const tangent_vector* i(directions.begin()); i != directions.end(); ++i) 

{ 

const point mirror_origin(*coordinateJine(*i, c)(arm_length)); 
const tangent_vector mirror_tangent(normalise(tangent_vector( 

mirror_origin, c, coordinate_direction))); 
const cached_worldline mirror(coordinate_line(mirror_tangent, c)); 

const geodesic outward_ray(connecting_null_geodesic(origin, mirror, 

armJength)— >second); 
const point reflection(*outward_ray(l.)); 

const pair<double, geodesio con(*connecting_null_geodesic(reflection, 

beam_splitter, two_arm_length)); 
const double&L affineJength(con. first); 
const geodesic&i inward_ray(con. second); 

plot(mirror, 0., affineJength); 

plot(outward_ray, 0., 1.); 
plot(inward_ray, 0., 1.); 

cout << "^^Mirror^direction^=J' << *(*i)[c] << " , Jime^experienced^ 
by,^beam,^splitter,^=J' << affineJength << "." << endl; 

cout << "^^^^Outward jay^tangent^=J' << *(*outward_ray.tangent(0.)) 
[c] << endl; 

cout << "^^^Jnward^ray^tangent^=J' << *(*inward_ray.tangent(0.))[c] 

<< endl; 
output << affineJength << endl; 

if (affineJength > max_affineJength) 
max_affineJength — affineJength; 

} 



plot(beam_splitter, 0., max_affineJength); 

} 



